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Abstract 

This  research  effort  analyzes  the  fundamental  dynamics  governing  a  satellite  with  a 
gravity  gradient  boom  and  a  tethered  balloon.  Satellites  that  use  gravity  gradient  booms  for 
passive  attitude  control  are  characterized  by  undamped  pitch  oscillations  and  no  roll  control.  The 
tethered  balloon  acts  as  a  high  drag  device  that  accounts  for  the  most  drag  on  the  satellite  system. 
By  attaching  a  drag  device,  the  system  resists  rolling  movements  while  also  damping  oscillations. 
This  could  potentially  be  a  cost  effective  method  for  increasing  satellite  stabilization.  The  goal  of 
this  research  is  to  model  the  dynamics  and  determine  the  feasibility  of  a  gravity  gradient 
stabilized  satellite  with  an  attached  balloon.  A  simulation  written  in  Matlab  analyzes  the  behavior 
of  such  a  satellite.  The  research  is  limited  to  circular  orbits  around  a  spherical  Earth  and  includes 
only  in-plane  motion  for  each  mass.  Stable  ranges  for  certain  tether  characteristics  are  found  for 
three  different  satellites. 
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A  NUMERICAL  ANALYSIS  OF  PASSIVE  ATTITUDE  STABILIZATION  USING 


A  TETHERED  BALLOON  ON  A  GRAVITY  GRADIENT  SATELLITE 

I.  Introduction 

This  research  effort  analyzes  the  fundamental  dynamics  governing  a  gravity  gradient  satel¬ 
lite  with  a  tethered  balloon  in  order  to  determine  the  feasibility  of  using  such  a  satellite 
concept  for  passive  attitude  stabilization.  Satellites  that  use  gravity  gradient  booms  for  pas¬ 
sive  attitude  control  are  characterized  by  undamped  oscillations  and  no  roll  control.  The 
tethered  balloon  acts  as  a  high  drag  device  that  accounts  for  the  most  drag  on  the  satel¬ 
lite  system,  which  causes  the  system  to  resist  rolling  movements  while  also  damping  pitch 
angle  oscillations.  This  could  potentially  be  a  cost  effective  method  for  increasing  satellite 
attitude  stabilization.  This  chapter  discusses  the  motivation  behind  the  analysis  of  such  a 
system.  Then,  the  objectives  of  this  research  are  listed,  the  approach  taken  to  achieve  those 
objectives  is  given,  and  the  limitations  of  this  research  are  discussed.  The  next  chapter  goes 
on  to  provide  a  more  adequate  background  for  this  satellite  concept.  Subsequent  chapters 
go  on  to  derive  the  equations  of  motion  for  the  satellite  concept,  detail  the  process  taken  in 
simulating  the  satellite  concept  in  orbit,  and  discuss  the  results  of  the  simulation  program. 

1.1  Motivation 

“Better,  Faster,  Cheaper”  is  a  common  catch-phrase  in  the  United  States  Air  Force 
(USAF).  In  a  service  centered  on  maintaining  technological  superiority  over  its  rivals,  this 
phrase  has  become  the  goal  for  many  Air  Force  engineers.  As  expressed  by  the  former 
Secretary  of  the  Air  Force  Dr.  James  G.  Roche,  “There  isn’t  enough  money  in  the  budget 
to  replace  everything  we  want  nor,  in  some  respects,  to  replace  everything  we  need.  We 
need  to  remain  committed  to  investing  wisely  in  our  future”  (7).  The  Air  Force  must  be 
cost  efficient  in  order  to  sustain  the  best-equipped  military  with  limited  resources. 

Furthermore,  the  actual  budget  for  the  Air  Force  space  program  is  quite  small  com¬ 
pared  to  the  Department  of  Defense  (DoD)  budget  and  even  the  Air  Force  budget.  The 
DoD  was  given  380  billion  dollars  in  2004,  from  which  93.5  billion  dollars  was  allocated  to 
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the  Air  Force  (5:4).  In  the  end,  the  Air  Force  allocated  about  1.8  billion  dollars  to  space- 
related  operations  and  31.2  billion  dollars  for  modernization,  a  portion  of  which  includes 
research,  development,  testing,  evaluation,  science,  technology,  and  procurement  of  space 
systems  (5:4).  This  is  a  high  estimate  for  space-related  programs  since  the  modernization 
budget  actually  includes  procurement  of  aircraft  and  armament.  While  space  technology  is 
a  high  priority,  Air  Force  space  programs  must  be  very  efficient  with  the  limited  funds  that 
it  receives. 

Commitment  to  advanced  technology  is  also  a  prime  characteristic  of  the  U.S.  Air 
Force.  As  stated  by  the  USAF  Future  Concepts  and  Transformation  Division: 

“The  purpose  of  Air  Force  innovation  is  to  rapidly  assess  and  implement  new 
ideas,  concepts,  and  technologies  to  field  the  best  capabilities  to  the  warfighter 
while  also  improving  the  associated  doctrine,  organization,  training,  materiel, 
leadership  and  education,  personnel,  and  facilities”  (14:21). 

The  Air  Force  was  founded  on  the  exploitation  of  new  technology  and  has  continued  to  be 
the  best  air  force  in  the  world  because  it  researches,  develops,  and  implements  the  newest 
technology. 

The  current  trend  in  satellites  is  a  shift  away  from  large,  multi-payload  satellites, 
which  carry  larger  production  and  launch  costs,  to  micro  and  nanosatellites  that  cost  much 
less  (8:806).  “Recently,  spacecraft  have  become  more  diverse,  with  the  largest  spacecraft 
now  complimented  by  new  systems  using  a  larger  number  of  smaller  spacecraft  in  low-Earth 
orbit  (LEO)”  (16:853).  Consequently,  the  smaller  size  and  lower  funds  drive  down  the  mass 
and  cost  budgets  of  the  satellite  subsystems.  To  keep  this  transition  effective,  the  satellite 
and  satellite  subsystems  must  maintain  a  high  degree  of  reliability.  In  addition,  increased 
subsystem  reliability  decreases  the  need  for  back-up  systems  which  in  turn  decreases  the 
total  satellite  cost  and  weight. 

Like  all  other  subsystems,  the  attitude  control  subsystem  must  meet  the  requirements 
determined  by  each  mission  payload  with  the  stricter  weight  and  cost  constraints  and  at 
the  same  degree  of  reliability.  In  fact,  the  attitude  control  subsystem  carries  additional 
importance  because  it  has  the  potential  for  decreasing  satellite  operation  costs  as  well  as 
production  costs.  With  greater  pointing  stability,  a  satellite  can  decrease  the  size  of  the 
communications  subsystem  by  decreasing  the  antenna  size  while  maintaining  a  low  level 
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of  bit  error.  With  a  totally  passive  attitude  control  system,  reliability  is  increased  and 
operations  costs  are  decreased.  A  passive  system  has  the  added  benefits  of  having  no 
fuel  requirements  and  very  little  power  requirements.  The  benefits  of  the  attitude  control 
subsystem  are  felt  by  all  satellite  subsystems,  decreasing  the  overall  cost  and  complexity  of 
the  satellite. 

Additionally,  few  low  budget  missions  have  the  financial  means  to  venture  beyond 
low-earth  orbit,  which  includes  up  to  500  kilometers  above  the  Earth.  Gravity  gradient 
and  aerodynamic  drag  torques  are  the  dominating  perturbations  at  these  altitudes.  Where 
many  satellites  must  counteract  the  effects  of  these  perturbations,  it  would  be  advantageous 
for  an  attitude  control  subsystem  to  use  these  torques  to  stabilize  the  satellite  and  meet 
the  mission  pointing  requirements.  A  satellite  using  passive  aerostabilization  decreases 
complexity,  weight,  and  cost  by  manipulating  these  perturbative  forces  rather  than  working 
against  them. 

1.2  Satellite  Concept 

A  satellite  concept  that  could  tie  together  all  of  these  issues  would  be  a  low  cost,  high 
performance,  and  high  reliability  passive  attitude  control  system  that  would  use  gravity 
gradient  and  drag  torques  in  low-Earth  orbit  to  produce  a  stable,  nadir-looking  attitude 
within  a  limited  time  frame.  The  proposed  answer  to  such  an  attitude  subsystem  would 
look  somewhat  like  Figure  1.1.  Essentially  starting  from  a  gravity  gradient  satellite,  a  main 
buss  mass  has  an  extended  gravity  gradient  boom  with  a  tip-mass  at  the  end.  Attached 
to  these  two  masses  is  a  high  drag  device,  potentially  a  rigidizable  balloon.  The  balloon  is 
attached  by  two  tethers  that  are  connected  to  the  tip-mass  and  main  buss  mass.  Figure  1.1 
shows  the  satellite  concept  in  an  initial  configuration,  but  it  is  expected  that  the  two  tethers, 
with  high  moduli  of  elasticity  and  elongation  percentages,  would  extend  the  balloon  mass 
further  behind  the  gravity  gradient  satellite.  A  deeper  explanation  of  this  configuration  is 
contained  in  Section  3.1. 
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1.3  Research  Objectives 

This  research  derives  and  tests  a  dynamical  model  for  a  gravity  gradient  satellite  that 
uses  a  tethered  balloon  for  completely  passive  attitude  stabilization.  The  purposes  of  this 
study  is  to  answer  the  following  questions. 

•  What  are  the  equations  of  motion  for  a  gravity  gradient  satellite  with  a  rigid  balloon 
attached  by  tethers? 

•  What  attitude  control  system  characteristics  can  be  altered  to  manipulate  the  effec¬ 
tiveness  and  efficiency  of  passive  gravity  gradient  attitude  control  with  aerostabiliza- 
tion? 

•  What  is  the  impact  of  certain  satellite  characteristics,  namely  modulus  of  elasticity 
and  damping  coefficient,  on  the  steady- state  attitude  of  this  satellite  concept? 

•  Where  are  the  regions  of  stability  for  these  satellite  characteristics? 

The  ideal  result  of  this  thesis  would  be  to  prove  or  disprove  the  feasibility  of  using 
a  gravity  gradient  boom  and  a  tethered  balloon  to  stabilize  the  attitude  of  a  satellite  by 
developing  and  testing  a  dynamical  model  for  the  satellite  system  through  a  computer 
program.  A  demonstration  of  a  stable  attitude  would  be  shown  through  the  damping  of  the 
gravity  gradient  pitch  angle  over  time.  This  is  the  basis  for  study  of  this  attitude  control 
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concept,  and  subsequent  design  studies  rely  on  whether  pointing  requirements  can  be  met 
at  the  base  research  level. 

1.4  Research  Approach 

In  order  to  accomplish  the  above  established  research  objectives,  this  study  takes  the 
following  steps: 

•  Review  all  relevent  information  surrounding  this  satellite  concept,  including  attitude 
control  techniques,  tethers,  aerostabilization,  and  any  other  pertinent  articles. 

•  Derive  the  equations  of  motion  for  a  gravity  gradient  satellite  with  a  tethered,  rigid 
balloon. 

•  Create  and  validate  a  computer  simulation  that  models  a  gravity  gradient  satellite 
with  aerostabilization. 

•  Characterize  the  impact  of  the  moduli  of  elasticity  and  damping  coefficients  of  the 
tethers  on  the  steady-state  attitude  of  a  satellite. 

•  Determine  the  regions  of  stability  for  these  satellite  characteristics. 

1.5  Research  Scope  and  Limitations 

The  research  objectives  allow  a  lot  of  room  for  simplifying  the  system  as  much  as 
possible  in  order  to  focus  on  the  dynamics  of  this  particular  concept;  however,  this  baseline 
study  must  be  thorough  enough  to  determine  whether  this  method  of  attitude  stabilization 
is  even  possible.  The  following  assumptions  limit  the  extent  of  knowledge  gaining  in  this 
research  effort  but  also  concentrate  all  efforts  to  produce  a  deeper  understanding  at  the  base 
level  of  research. 

•  Assume  two-body  equations  of  motion  for  the  satellite  system  center-of-mass  (COM) 

•  Assume  circular,  low-Earth  orbit  about  a  spherical  Earth 

•  Assume  gravity,  drag,  tether-spring,  and  tether-damper  are  the  only  forces  acting  on 
the  tip-mass,  main  satellite  bus,  and  tethered  balloon 
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•  Assume  motions  of  the  tip-mass,  main  satellite  bus,  and  tethered  balloon  are  limited 
to  orbital  plane 

•  Assume  rigid  masses  and  point  masses  for  tip-mass,  main  satellite  bus,  and  tethered 
balloon 

•  Assume  tethers  are  not  compressed  and  are  massless 

•  Assume  rigid  and  massless  gravity  gradient  boom 

•  Assume  uniform  density  for  atmosphere  at  a  specific  altitude 

•  Focus  on  motion  of  the  satellite  after  deployment  of  the  gravity  gradient  boom  and 
tethered  balloon 

These  assumptions,  while  limiting  the  scope,  greatly  simplify  the  dynamics  of  the  satellite 
system.  In  most  cases,  the  calculations  and  equations  involved  with  this  study  are  less 
complicated  and  less  prone  to  error.  Further  discussion  of  these  assumptions  and  issues  is 
contained  in  Section  3.1. 
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II.  Literature  Review 


The  Literature  Review  Chapter  discusses  pertinent  background  material  that  is  helpful  in 
understanding  this  satellite  concept.  Since  the  satellite  configuration  that  is  being  studied 
involves  a  gravity  gradient  boom  and  tethered  balloon  used  for  passive  attitude  control,  this 
chapter  will  describe  attitude  control  techniques,  gravity  gradient  booms,  tethers  in  space, 
and  rigidizable  balloons.  Past  research  involving  passive  aerostabilization  is  also  discussed 
since  it  is  the  closest  configuration  to  the  presently  studied  satellite  concept.  While  this  is 
not  a  comprehensive  description  of  all  satellite  topics  necessary  to  understand  this  study, 
this  chapter  attempts  to  explain  and  discuss  all  immediately  relevant  information. 

2.1  Attitude  Control  Techniques 

As  defined  by  Space  Mission  Analysis  and  Design  (SMAD),  “the  attitude  determi¬ 
nation  and  control  subsystem  measures  and  controls  the  spacecraft’s  angular  orientation 
(pointing  direction)  or  its  orientation  and  linear  velocity”  (16:302).  Pointing  control  is  nec¬ 
essary  for  several  different  reasons  and  at  different  performance  specifications,  depending 
on  the  mission  of  the  satellite.  Reasons  for  attitude  control  include  orbit  insertion,  initial 
attitude  stabilization,  normal  station-keeping,  and  special  or  contingency  slew  maneuvers 
(16:356).  These  different  operation  modes  can  result  in  totally  different  attitude  control  sys¬ 
tems.  Additionally  under  the  normal  station-keeping  mode,  satellites  are  subject  to  cyclic 
and  secular  disturbance  torques.  Cyclic  disturbance  torques  happen  on  a  periodic  basis, 
while  secular  torques  have  linear  effects,  like  the  effect  of  solar  pressure.  Furthermore,  grav¬ 
ity  and  drag  dominate  the  disturbance  torques  at  LEO  while  solar  pressure  and  third-body 
effects  have  greater  influence  in  high-Earth  orbit  (HEO).  Determining  the  primary  distur¬ 
bance  torques  in  the  mission  orbit  can  help  the  design  engineers  focus  on  certain  issues.  All 
of  these  requirements  and  issues  drive  the  design  of  the  attitude  control  system. 

In  designing  the  attitude  control  system,  techniques  come  in  the  form  of  passive  sys¬ 
tems,  spin-control  systems,  and  three-axis  control  systems.  A  gravity  gradient  boom  is 
an  example  of  a  passive  system  and  is  described  further  in  Section  2.2.  Passive  systems 
normally  have  decreased  accuracy  but  are  less  expensive  and  complex.  “They  consume  no 
power,  require  no  hardware  for  sensing  or  actuation,  and  make  no  demands  on  software” 
(6:282).  The  opposite  is  true  of  active  systems.  Thrusters  and  magnetic  torquers  are  ex- 


2-1 


amples  of  active  systems.  Active  systems  use  either  of  two  concepts,  zero  momentum  or 
momentum  bias.  In  zero-momentum  systems,  “reaction  wheels  respond  to  disturbances  on 
the  vehicle”  (16:362).  Momentum  bias  systems  use  one  spinning  wheel  to  give  the  satel¬ 
lite  gyroscopic  stiffness  and  to  control  the  attitude  by  torquing  the  wheel.  Spin  stabilized 
satellites  passively  resist  disturbance  torques  by  using  gyroscopic  stability,  but  extra  fuel  is 
needed  to  reorient  the  satellite  once  it  begins  spinning. 

In  determining  which  attitude  control  technique  to  use  on  a  satellite,  it  is  important 
to  look  at  the  required  performance  specifications  in  terms  of  accuracy,  range,  jitter,  drift, 
and  settling  time  (16:357).  For  example,  “In  many  modern  applications,  the  attitude  errors 
permitted  are  so  small  that  only  a  fully  automatic  (fully  ‘active’)  control  system  can  meet 
mission  specifications”  (6:281).  These  requirements  may  come  from  different  subsystems, 
for  example  the  payload  or  the  communications  subsystem.  Accuracy  and  settling  time  are 
the  focus  of  this  feasibility  study  since  this  satellite  concept  is  for  low-budget  microsatellites. 
More  stringent  requirements  would  be  used  for  further  studies  and  designs. 

2. 2  Gravity  Gradient  Satellites 

The  prime  example  for  passive  attitude  control  is  gravity  gradient  stabilization.  Grav¬ 
ity  gradient  satellites  exploit  “some  naturally  occurring  force  field  [Earth’s  gravity  field]  to 
provide  the  desired  stabilizing  torque”  (6:282).  Figure  2.1  shows  an  example  of  a  satellite 
that  uses  gravity  gradient  stabilization. 

In  designing  a  gravity  gradient  stabilized  satellite,  it  is  important  for  the  moments 
of  inertia  to  follow  certain  criteria.  First,  the  inertia  ratio  parameters  k\  and  k%  must  be 
defined.  These  ratios  essentially  make  the  design  a  function  of  two  variables,  rather  than 
three.  In  these  equations,  I\  is  the  moment  of  inertia  aligned  with  the  satellite’s  velocity 
vector,  I2  is  normal  to  the  orbital  plane,  and  1%  is  aligned  with  the  vertical  axis. 

*1  =  ^  (2.1) 

h 

k3  =  (2.2) 
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Figure  2.1  GEOS-II  Spacecraft  (6:328) 

The  regions  of  stability  can  then  be  seen  in  Figure  2.2.  Hughes  outlines  the  process  of 
determining  the  stability  regions  for  gravity  gradient  satellites  in  chapter  nine  of  his  book, 
Spacecraft  Attitude  Dynamics  (6).  From  this  figure,  it  can  be  seen  that  only  certain  satellite 
configurations  are  acceptable  for  gravity  gradient  stabilization.  The  upper-right  region  is 
populated  by  satellites  that  have  the  minor  axis  in  the  vertical  direction,  the  intermediate 
axis  pointing  in  the  velocity  direction,  and  the  major  axis  pointing  normal  to  the  orbit.  The 
lower  region,  called  the  DeBra-Delp  region,  shows  the  stability  of  rigid  bodies  whose  minor 
axis  is  normal  to  the  orbit  and  whose  major  axis  is  tangent  to  the  orbit  (6:301).  The  satellite 
concept  under  study  uses  the  region  of  stability  described  in  the  upper  right  of  Figure  2.2. 

Although  these  stable  regions  exist  in  theory,  there  are  two  problems  with  gravity 
gradient  stabilization.  First,  gravity  gradient  torques  decrease  with  increased  altitude.  This 
makes  it  necessary  to  use  inertia  augmentation,  which  means  increasing  the  moments  of 
inertia  as  much  as  possible.  This  can  be  done  by  placing  long,  slender  rods  along  the 
principal  axes.  Effectively,  moments  of  inertia  are  increased,  which  also  increases  gravity 
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Liapunov  Stable  and  Statically  Stable 


Statically  Unstable 
igogogp1:  but  Gyrically  Stabilized 


Figure  2.2  Stability  Regions  for  Gravity  Gradient  Satellites  in  Circular  Orbit  (6:297) 


gradient  restoring  torques.  This  concept  is  depicted  in  Figure  2.3.  The  second  issue  in 
dealing  with  gravity  gradient  stabilization  is  effective  damping  or  energy  dissipation.  In 
order  to  dampen  out  the  librations,  or  oscillations  in  the  pitch  direction,  the  system  must 
dissipate  energy.  Options  for  energy  dissipation  include  magnetic-hysteresis  rods,  spherical 
tip  dampers,  and  boom  articulation.  This  damping  is  also  a  major  focus  for  the  satellite 
concept  currently  under  investigation. 


2.3  Tethers  in  Space 

The  idea  of  tethers  in  space  has  been  around  since  1895  (3:7).  Since  then,  the  range 
of  applications  for  space  tethers  has  been  ever-growing.  The  usefulness  of  tethers  span 
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Figure  2.3  Inertia  Augmentation  (6:315) 


all  elements  of  the  satellite  system  architecture.  For  example,  tethered  probes  can  take 
samples  and  measurements  at  different  altitudes  without  having  to  maneuver  the  main 
satellite  or  spacecraft.  This  application  is  shown  in  Figure  2.4.  This  of  course,  eliminates 


Figure  2.4  Tethered  Probe  Taking  Measurements  at  Different  Altitude  than  Orbiter  (3:26) 

the  need  for  extra  maneuvering  fuel  and  expends  the  field  of  measurement  around  the  main 
satellite.  Tethers  also  have  the  potential  to  collect  electric  power  through  a  metallic  tether 
when  travelling  through  the  Earth’s  geomagnetic  field.  The  same  essential  concept  makes 
it  possible  to  produce  thrust  when  electricity  is  applied  to  a  metallic  tether.  These  two 
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concepts  are  depicted  in  Figure  2.5.  Tethers  have  even  more  applications  involving  only 


Figure  2.5  Metallic  Tethers  Travelling  through  Earth’s  Geomagnetic  Field  (a)Producing 
Thrust  (b)Producing  Electric  Power  (3:25) 

space  transportation.  Concepts  like  the  space  tether  elevator,  depicted  in  Figure  2.6,  could 
be  used  on  the  Earth  as  well  as  the  Moon.  Tether  elevators  attach  to  either  the  Earth  or 


Figure  2.6  Moon  Based  Space  Elevator  (3:30) 

Moon  and  extend  into  normal  orbital  altitudes.  Satellites  using  this  system  would  travel  up 
the  tether  elevator  and  would  be  implanted  into  orbit  using  much  less  fuel  than  a  typical 
ground- launched  rocket.  Space  escalators,  shown  in  Figure  2.7,  use  similar  concepts  to  space 
elevators  but  applied  to  several  tether  satellites  orbiting  at  different  altitudes.  Satellites 
attach  at  one  end  of  the  first  tether  satellite,  travel  to  the  opposite  end  of  the  tether,  release 
from  that  tether  satellite,  attach  to  the  next  tether  satellite  orbiting  further  away,  and 
repeat  that  cycle  until  the  desired  orbit  is  reached.  Since  the  tether  systems  effectively 
transfer  energy  to  the  payload  satellites  to  boost  it  to  higher  altitudes,  the  orbiting  tethers 
need  some  kind  of  station- keeping  subsystem  to  keep  each  tether  in  its  respective  orbit. 

There  are  many  hindrances  to  the  implementation  of  the  above-mentioned  tether  sys¬ 
tems.  Strength,  for  one,  is  a  material  specification  that  limits  the  length  of  the  tethers 
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Figure  2.7  Space  Escalator  (3:29) 


required  for  these  applications.  Specifically,  there  are  no  known  materials  that  can  be  pro¬ 
duced  to  the  lengths  necessary  for  a  space  elevator  (at  least  35,800  kilometers  to  reach 
geostationary  orbit).  A  tether  at  this  length  could  not  withstand  its  own  weight  and  def¬ 
initely  not  the  added  weight  of  a  satellite.  This  mainly  impacts  tethers  extending  from 
Earth  since  the  gravitational  force  acting  on  a  tether  at  even  the  lowest  Earth  orbits  is  large 
enough  to  “exceed  the  break  lengths  on  Earth  by  an  order  of  magnitude”  (3:36).  Mass  is 
another  issue  as  many  applications  of  space  tethers  require  lengths  in  units  of  kilometers. 
This  obstacle  is  crossed  by  using  several  braided  fibers  rather  than  one  long  strand.  This 
method  increases  strength  with  lower  densities.  Debris  and  micrometeorites  pose  another 
threat  since  a  great  majority  of  which  have  diameters  comparable  to  the  diameters  of  space 
tethers  and  travel  at  an  average  relative  speed  of  20  kilometers  per  second  (3:39).  While 
the  threat  is  small  for  tethers  smaller  than  100  kilometers,  tether  tapes  have  been  proposed 
for  decreasing  the  dangerous  impact  from  small  space  debris. 

2.4  Rigidizable  Balloons 

Rigidizable,  inflatable  balloon  systems  have  been  studied  since  the  beginning  of  the 
National  Aeronautics  and  Space  Administration  (NASA),  but  research  is  beginning  to  in¬ 
crease  due  to  modern  missions  which  require  even  larger  on-orbit  systems,  like  large  aperture 
optical  telescopes,  solar  arrays,  and  aerobrakes.  Rigidizable  balloons  are  intended  to  replace 
typical,  rigid  mechanical  structures.  The  balloons  are  flexible  when  not  inflated  but  become 
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stiff  when  inflated  to  a  predesigned  size.  Such  systems  are  packaged  uninflated  prior  to 
launch  in  order  to  fit  what  would  be  systems  larger  than  the  launch  vehicle  envelope.  After 
the  satellite  is  launched,  the  balloon  inflates  to  a  predesigned  shape.  Rigidizable  balloons 
derive  their  motivation,  in  large  part,  from  the  same  motivation  for  this  research  concept  - 
smaller,  more  light-weight  systems.  The  cost  of  a  satellite  increases  as  its  size  and  weight 
increases,  but  inflatable  balloons  can  be  fifty  percent  lighter  and  twenty- five  percent  smaller 
than  standard  mechanical  systems  (12:2-3). 

The  system  has  high  reliability  because  it  does  not  require  complex,  mechanical  com¬ 
ponents  or  joints.  Hinges  and  joints  on  mechanical  systems  have  the  potential  to  cold  weld, 
making  it  impossible  to  extend  or  deploy  a  part  of  the  satellite.  On  the  other  hand,  in¬ 
flatable  balloons  continue  to  increase  the  pressure  and  force  of  deployment  in  the  event  the 
system  somehow  resists  deployment  (12:2-3).  The  only  failure  point  on  an  inflatable  system 
is  the  trigger  to  initiate  inflation.  Also,  inflatable  systems  can  be  predesigned  to  inflate  to 
any  shape  and  can  even  be  packaged  to  almost  any  shape.  This  flexibility  can  be  crucial 
when  trying  to  fit  a  satellite  system  within  a  launch  vehicle’s  envelope. 

Slightly  different  from  rigidizable,  inflatable  structures  are  inflatable  structures,  which 
require  additional  gas  in  order  to  maintain  the  structural  shape  in  the  event  of  leaks.  In¬ 
flatable  structures  have  been  the  focus  of  many  research  efforts  even  though  rigidizable, 
inflatable  structures  do  not  require  extra  gas  and  do  not  cause  attitude  torques  in  the  case 
of  punctures  in  the  balloon,  possibly  from  micrometeorite  impacts.  Any  leaks  in  inflatable 
devices  could  cause  unwanted  thrust  due  to  gas  being  expelled  through  small  holes  in  the 
structure.  The  inflatable  structure  also  adds  weight  to  the  total  satellite  because  of  the  need 
of  additional  gas,  gas  tanks,  and  regulators. 

An  example  of  an  inflatable,  rigidizable  structure  is  the  L’Garde  Inflatable  Space  Truss, 
shown  in  Figure  2.8  .  The  truss  is  made  out  of  a  “composite  material  that  are  impregnated 
with  a  water  soluble  resin”  (12:2-11).  The  water  in  the  composite  evaporates  in  the  vacuum 
of  space,  stiffening  the  resin  material  and  the  structure  as  a  whole.  This  system  had  a  final 
weight  of  1.917  kilograms  and  a  length  of  60.1  inches,  although  the  packing  volume  was  only 
1953  cubic  inches  (12:2-12). 
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Figure  2.8  L’ Garde  Inflatable  Space  Truss  (12:2-12) 


As  inflatable,  rigidizable  systems  are  still  under  development,  there  are  several  prob¬ 
lems  that  this  technology  must  overcome.  First,  it  is  difficult  to  test  the  effects  of  space  on 
these  systems  while  they  are  on  the  ground.  For  example,  the  effect  of  gravity  overshadows 
the  force  seen  in  space.  Since  these  forces  cannot  be  properly  tested,  the  guidance  and 
control  of  the  satellite  gains  additional  demands.  Gravity  also  plagues  the  fabrication  of 
the  material  and  the  structure  on  the  ground.  The  effects  of  gravity  can  make  it  difficult  to 
achieve  uniformity  in  the  material.  This  can  cause  a  distortion  in  the  shape  when  deployed 
in  space.  These  problems  will  hopefully  be  overcome  by  creating  modeling  and  simulation 
software  that  can  more  accurately  test  these  structures  and  materials  in  a  simulated  space 
environment  without  actually  being  in  space. 

2.5  Past  Research 

NASA  was  the  first  to  develop  and  demonstrate  a  totally  passive,  aerodynamically 
stabilized  satellite  for  LEO  in  the  years  1993  to  1997.  The  NASA  Goddard  Space  Flight 
Center  developed  the  Passive  Aerodynamically  Stabilized  Magnetically  Damped  Satellite 
(PAMS),  shown  in  Figure  2.9,  as  the  first  flight  experiment  to  demonstrate  the  concept  of 
passive  aerostabilization.  This  technology  demonstration  was  required  so  that  this  concept 
could  be  used  on  another  satellite,  the  Gravity  and  Magnetic  Earth  Surveyor  (GAMES), 
which  would  be  flown  in  order  to  evaluate  the  gravity  field  of  the  Earth  to  a  high  precision. 
Passive  aerostabilization  was  seen  as  a  “low-cost,  low- weight,  long-lifetime  option”  for  the 
satellite  (9:228). 
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Figure  2.9  PAMS  Schematic  (9:228) 

Figure  12.  7  Day  Cone  Angle  History  -  Modified  ICs 


FREEMOL  Simulated  Rs&ulEs 


Elapsed  Qrbiss  rrom  Deploymofil 

Figure  2.10  PAMS  Results  (15) 

PAMS  used  a  configuration  similar  to  a  shuttlecock  where  the  center  of  mass  was 
placed  forward,  in  the  ram  direction,  of  the  center  of  pressure,  while  magnetic  hysteresis 
rods  were  used  for  rate  damping.  Aerodynamic  stabilization  was  chosen  as  the  primary  sta¬ 
bilization  force  because  as  Psiaki  explains,  “At  altitudes  below  400  km,  aerodynamic  drag 
torque  tends  to  overwhelm  the  gravity  gradient  torque  for  practical  lightweight  deployable 
boom  designs”  (13:2).  NASA’s  method  for  attitude  stabilization  avoided  the  gravity  gradi¬ 
ent  disturbance  torque  by  designing  the  satellite  to  have  equal  moments  of  inertia.  NASA 
produced  the  computer  program  Free  Molecular  Aerodynamic  Satellite  Attitude  Dynamic 
Simulation  (FREEMOL)  to  predict  and  numerically  confirm  the  feasibility  of  aerostabi- 
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lization.  This  simulation,  shown  in  Figure  2.10,  predicted  worst-case  angular  rates  of  0.1 
degress  per  second  in  yaw  and  pitch  and  0.05  degrees  per  second  in  roll. 


Figure  2.11  NASA  Flight  of  PAMS  (15) 


PAMS  was  flown  on  the  Space  Shuttle  Mission  STS-77  on  May  1996.  The  PAMS 
flight  measurements  met  complications  that  prevented  validation  and  calibration  of  the 
FREEMOL  software,  but  the  cone  angle  estimates  made  by  space  shuttle  astronauts  during 
rendezvous  verified  that  the  satellite’s  attitude  stabilized.  During  flight,  the  instrument 
intended  to  measure  the  satellite’s  attitude  failed.  Although  this  portion  of  the  PAMS 
flight  was  a  failure,  the  overall  project  was  still  deemed  a  success.  While  PAMS  is  a  similar 
passively  stabilized  system  developed  and  tested  by  NASA,  the  satellite  concept  currently 
under  study  has  never  been  tried. 
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III.  Dynamics 

This  section  will  outline  and  describe  the  fundamental  mathematics  that  are  used  to  model 
the  dynamics  of  this  satellite  concept.  The  chapter  first  describes  the  satellite  concept  and 
the  assumptions  made.  Next,  the  chapter  defines  the  coordinate  frames  used  in  this  study. 
Finally,  the  two  body  equation  of  motion  is  derived  as  well  as  the  equations  of  motion  for 
this  satellite  concept. 

3.1  Basic  Configuration  and  Assumptions 

The  distance  from  the  gravity  gradient  tip-mass  to  the  center  of  mass  is  defined  as 
’b’  whereas  the  distance  to  the  main  satellite  bus  from  the  tip-mass  is  defined  as  ’a’.  The 
length  of  each  tether  when  taut  is  defined  as  d\  for  tether  one  and  for  tether  two.  The 
moduli  of  elasticity  for  each  tether  are  '>E\  and  ’’EC  and  the  damping  coefficients  are  ’ci’ 
and  ’c2?,  respectively.  Also,  the  tip-mass  is  defined  as  mass  one  or  mi,  the  main  bus  mass 
is  m2,  and  the  balloon  mass  is  m3.  Tether  one  connects  the  tip-mass  to  the  balloon  mass, 
and  tether  two  connects  the  main  bus  mass  to  the  balloon  mass.  These  characteristics  are 
depicted  in  Figure  3.1. 

For  this  study,  the  tension  in  the  tether  is  mathematically  modelled  as  a  simple  spring 
with  a  modulus  of  elasticity.  Normally,  moduli  of  elasticity  are  in  units  of  force  per  cross- 
sectional  area  of  the  spring.  In  this  case,  the  moduli  of  elasticity  have  units  of  Newtons  per 
millimeter  squared,  but  the  cross-sectional  area  is  assumed  to  be  one  millimeter  squared. 
In  order  to  make  use  of  this  assumption,  the  moduli  of  elasticity  are  in  units  of  Newtons 
in  the  equations  of  motion  and  in  the  simulation  program.  However,  this  characteristic 
is  listed  in  force  per  unit  area  in  the  rest  of  this  paper.  The  damping  coefficient  serves 
as  a  model  for  energy  dissipation  in  the  spring.  Energy  dissipation  is  the  factor  that  will 
cause  the  librations,  or  oscillations  due  to  gravity  gradient  torque,  to  dampen.  The  damping 
coefficients  are  in  units  of  Newton  seconds  per  meter.  This  is  regardless  of  the  cross-sectional 
area,  and  these  are  the  units  used  for  this  characteristic  in  the  dynamical  model. 

The  center  of  mass  of  the  satellite  system  is  approximated  as  residing  on  the  gravity 
gradient  boom  at  all  times.  The  gravity  gradient  boom  and  tethers  are  assumed  to  have 
zero  mass.  The  tethered  balloon  is  assumed  to  have  a  mass  much  smaller  than  the  tip-mass. 


3-1 


Figure  3.1  Depiction  of  Satellite  Concept 

By  definition,  the  satellite’s  main  bus  has  much  more  mass  than  the  tip-mass.  Also,  the 
satellite  is  modelled  as  if  the  gravity  gradient  boom  were  already  deployed.  In  addition, 
the  tethered  balloon  is  assumed  to  have  been  deployed  and  inflated.  The  dynamics  in  this 
model  begin  as  if  the  satellite  system  depicted  in  Figure  3.1  is  deployed  as  shown.  As  the 
steady  state  stabilization  is  the  focus  of  this  baseline  study,  it  is  not  necessary  to  consider 
the  system’s  response  in  the  transient  mode.  In  fact,  the  mechanism  for  deploying  both  the 
gravity  gradient  boom  and  tethered  balloon  is  not  within  the  scope  of  this  study.  These 
assumptions  are  in  addition  to  those  listed  in  Section  1.5. 

3.2  Coordinate  Frames 

To  keep  the  model’s  dynamics  as  simple  as  possible,  this  study  uses  a  local  coordinate 
frame  to  the  satellite  as  shown  in  Figure  3.1  of  Section  3.1.  This  non-inertial  coordinate 
frame’s  first  two  axes  are  in  the  orbital  plane  with  the  origin  at  the  center  of  mass.  The  first 
axis,  points  in  the  direction  of  the  satellite’s  velocity  vector,  the  second  axis,  5y ,  points 
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in  the  radial  direction  away  from  Earth,  and  the  third  axis,  Sz,  points  out  of  the  orbital 
plane,  completing  the  right-handed  coordinate  frame. 

The  radius,  rc,  of  the  satellite’s  orbit  is  measured  from  the  center  of  the  Earth  to 
the  center  of  mass  of  the  satellite  system.  For  this  study,  the  altitude  is  assumed  to  be 
200  kilometers  above  the  Earth’s  surface  ( REarth  —  6375 km).  This  project  also  assumes  a 
circular  orbit  about  a  spherical  Earth. 

k,  north,  u 


J 

l.T 


Figure  3.2  Geocentric  Equatorial  Coordinate  Frame 

The  Geocentric  Equatorial  (GCE)  coordinate  frame,  which  is  inertially  fixed  in  space 
and  shown  in  Figure  3.2,  is  used  in  this  study  in  defining  the  two-body  equations  of  motion 
in  Section  3.3.  The  first  axis  points  toward  the  vernal  equinox,  the  second  axis  is  normal  to 
the  first  and  in  the  equatorial  plane,  and  the  third  axis  is  normal  to  the  first  two  axes. 

3.3  Two-Body  Equation  of  Motion 

The  two-body  equation  of  motion  mathematically  describes  the  motion  of  a  satellite. 
The  two  bodies  in  this  system  are  the  Earth,  with  mass  ra#,  and  a  satellite,  with  mass 
ras,  both  assumed  to  be  point  masses  as  shown  in  Figure  3.3.  With  respect  to  the  inertial 
frame,  the  positions  of  each  mass  are  Re  and  Rs  while  the  accelerations  are  Re  and  Rs. 
The  position  of  the  satellite  with  respect  to  the  Earth  is  defined  as  f*0,  where: 

r0  =  Rs  —  Re  (3.1) 
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Figure  3.3  Inertial  Frame  for  Two-Body  Equations  of  Motion 

The  acceleration  of  the  satellite  with  respect  to  the  Earth  is  determined  by  taking  the  second 
derivative  of  the  position  vector,  f*Q. 


r0  =  Rs  —  Re 


(3.2) 


Using  Newton’s  second  law,  F  —  ma,  the  only  force  acting  on  each  mass  is  the  gravitational 
force  due  to  the  other  mass  and  is  defined  by  the  following  equations. 


Fs  =  - 


FE  =  - 


Gms  mE 


Rs  —  Re 


Gms  mE 


(Rs  Re) 


[RE  —  Rs) 


(3.3) 


(3.4) 


|  Re  ~  Rs\ 

In  Equations  3.3  and  3.4,  G  is  defined  as  the  gravitation  constant.  After  substituting  the 
above  equations  into  Newton’s  second  law  and  replacing  a  with  RE  and  Rs ,  Equations  3.5 
and  3.6  are  found. 

(R,  -  Re)  (3.5) 


ms  iL  =  — 


Gms  mE 


mE  Re  = 


!  >  —  Re 
Gms  mE 


Re  —  Rs 


(■ Re 


-R, 


(3.6) 
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Substituting  these  into  Equation  3.2,  the  equation  of  motion  is  found  to  be: 


Tn  =  -- 


G  (ms  +  mE) 


Rs  —  Re 


(^Rs  -  Re^J 


(3.7) 


To  further  simplify  this  equation  of  motion,  the  standard  gravitational  parameter  is  defined 
in  Equation  3.8. 

fi  =  G  (ms  +  mE)  (3.8) 


This  quantity  is  known  more  accurately  than  the  separate  constituent  values.  For  man¬ 
made  satellites,  (i  is  approximated  as  /i  =  Gttie  since  ms  is  much  smaller  than  m#.  After 
substituting  /i  and  Equation  3.1  into  Equation  3.7,  the  final  two-body  equation  of  motion 
is: 

ro  =  --^Gr0  (3.9) 

\r0\ 

The  Geocentric  Equatorial  (GCE)  coordinate  frame,  described  in  Section  3.2,  and  the 
two-body  equation  of  motion  are  used  to  propagate  the  position  and  velocity  vectors  of  the 
satellite’s  center  of  mass  forward  in  time.  The  derivation  of  the  two-body  equation  of  motion 
is  fundamentally  extracted  from  and  further  explained  in  Section  2.2  of  Spacecraft  Dynamics 
by  Dr.  William  Wiesel  (17). 


3.4  Satellite  Concept  Equations  of  Motion 

The  two-dimensional  configuration  for  the  studied  satellite  concept  is  the  most  straight¬ 
forward  case  on  which  to  determine  the  feasibility  of  this  system.  This  case  is  based  on  a 
circular  orbit  defined  by  two-body  motion  and  includes  the  force  of  drag.  The  local  coordi¬ 
nate  frame  used  in  this  case  is  described  in  Section  3.2.  The  process  used  to  determine  the 
satellite  concept’s  equations  of  motion  is  known  as  Lagrange’s  Equations  of  Motion.  This 
process  was  selected  for  its  simplicity  and  ability  to  factor  in  additional  non-conservative 
forces  without  altering  the  basic  equations  of  motion,  which  only  include  conservative  forces. 

The  equations  of  motion  for  this  satellite  concept  use  the  generalized  coordinates  fe, 
5y ,  and  6.  These  are  shown  in  Figure  3.1  of  Section  3.1.  There  are  only  three  coordinates 
because  there  are  three  masses  with  a  total  of  only  three  degrees  of  freedom.  The  number  of 
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degrees  of  freedom  is  found  by  subtracting  the  number  of  constraints  from  three  times  the 
number  of  masses.  From  the  assumptions,  none  of  the  masses  move  out  of  the  orbital  plane, 
which  equates  to  one  constraint  for  each  mass  or  three  total  constraints.  The  tip-mass  and 
main  bus  mass  can  move  in  the  Sx  and  Sy  directions,  but  both  are  fixed  in  distance  from 
the  center  of  mass  which  accounts  for  two  more  constraints  for  the  total  system.  Also,  the 
tip-mass  and  the  main  satellite  bus  are  fixed  rigidly  to  each  other,  accounting  for  the  sixth 
constraint.  Therefore,  three  masses  start  with  nine  total  degrees  of  freedom  from  which  six 
constraints  are  subtracted. 


The  balloon  is  essentially  free  to  move  in  the  5x  and  5y  directions  within  the  bounds 
of  the  two  tethers  so  Sx  and  Sy  are  a  good  choice  for  the  first  two  generalized  coordinates. 
The  tip-mass  and  main  satellite  bus  are  fixed  to  each  other  and  pivot  about  the  center  of 
mass.  This  is  the  same  as  a  pendulum  so  an  angle,  0,  is  a  typical  choice  for  the  generalized 
coordinate.  The  coordinate  Sy  is  measured  from  the  center  of  mass  of  the  satellite  system 
radially  away  from  the  center  of  the  Earth.  The  coordinate  Sx  is  measured  in  the  ram 
direction  emanating  from  the  center  of  mass  of  the  satellite  system.  The  pitch  angle,  0,  is 
measured  from  the  Sy  axis  to  the  gravity  gradient  boom,  b.  This  angle  is  defined  to  be 
positive  as  rotated  clockwise  about  the  negative  Sz  axis. 


To  begin  the  process  in  determining  the  equations  of  motion  for  this  system,  the 
position  of  each  mass  is  determined  in  terms  of  the  local  reference,  or  body,  frame.  Since 
all  masses  are  in  the  orbital  plane,  the  position  and  velocity  components  along  the  third 
axis  are  zero.  The  three  masses  (tip-mass,  main  satellite  bus,  and  balloon)  have  respective 
position  vectors: 


Sri  = 

Sr2  = 


'  Sx  i ' 
Jvi. 
'  Sx  2' 
,Sy2, 


'  b  sin(0) 

,  rQ  +  b  cos (0) , 
—a  sin(0) 
Kr0  —  a  cos (0) 


(3.10) 

(3.11) 


Sr 


3 


*  Sx  3 
V  <%3 


^  Sx  \ 

^rQ  +  Sy  J 


(3.12) 
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The  velocity  vectors  in  terms  of  the  body  frame  are  determined  by  taking  the  first  derivative 
of  the  position  vectors  with  respect  to  the  inertial  GCE  frame  according  to  Equation  3.13. 


I 

TB 


-4  B 
TB 


+  (w  x  rB) 


(3.13) 


Using  Equation  3.13  with  uj  equal  to  the  angular  velocity  of  the  satellite  system  about  the 
Earth, 

(3.14) 

the  velocities  of  each  mass  are  calculated  to  be: 


Sr  i 


Sr2  = 


/bO  cos (9)  —  uo  (ro  +  b  cos (9))  \ 
\  —  b6  sin(0)  +  uo  b  sin(0)  / 

'  —a  9  cos (6)  —  u  (ro  —  a  cos (0)) ' 


6r3  = 


(3.15) 


(3.16) 

(3.17) 


a  6  sin(0)  —  uo  a  sin (6)  J 

'  Sx-u;  (r0  +  Sy)  \ 

5y  +  uo  Sx  ) 

Next,  the  kinetic  energy  of  the  system  is  found  by  summing  the  kinetic  energies  of  each 
mass. 

T  =  ^  mi  ^5ri  •  5r^  +  ^  m2  fsr2  •  Sr2 )  +  ^  m3  fsr3  •  Sr3)  (3.18) 


The  potential  energy  of  the  system  is  found  in  similar  fashion  by  summing  the  gravitational 
potentials  acting  on  each  mass. 


P  rrii  i 

um2  i 

um3 

Sri 

8t2 

Sr  3 

(3.19) 


The  denominators  of  the  potential  energies  in  Equation  3.19  are  approximated  using  the 
binomial  theorem  expanded  to  order  three. 

(a  +  /3)n  —  an  +  nan~l  (3  +  —  n  (n  —  1)  an~ 2  /32  +  ...  (3.20) 
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By  making  the  following  substitutions: 


OL  =  Tn 


(3.21) 


P  =  (2  ray  +  y2  +  a;2) 


(3.22) 


where  x  and  y  are  place  holders  for  the  corresponding  components  of  the  positions  of  each 
mass,  Sri ,  this  approximation  is  carried  out  with  the  following  results: 


Sr  1 


-1 


+ 


'  o 

3 


=  rn  1 - ^7  f2  r0  b  cos (6)  +  62  cos(0)2  +  b 2  sin(#)2 

2  ry5  V 


8r05 


^2  rQ  b  cos (0)  +  62  cos(#)2  +  62  sin(#)2^  +  7?  (3) 


(3.23) 


-1 


—  rQ  1 - f — 2  rQ  a  cos (0)  +  a2  cos(0)2  +  a2  sin(0)2) 

2  Tvy  '  ' 


+  ... 


+ 


8r05 


(—2  rQ  a  cos (9)  +  a2  cos(0)2  +  a2  sin(0)2^  +  7?  (3) 


(3.24) 


Sr3 


-1 


-r  -1  _ 
—  1  n 


7^7  (2  ra  5y  +  5y2  +  5x2^j  +  (2  rD  5y  +  5y2  +  &c2)  +  •&  (3)  (3.25) 


The  results  of  Equations  3.18  and  3.19  are  used  to  define  the  Lagrangian. 


L  —  T  —  V 


(3.26) 


The  Lagrangian  equations  of  motion  are  found  for  n-generalized  coordinates  by  Equation 
3.27. 

d_  (dL\  _  f  d  L  \  _ 


Qki  k  lj  2,  ...77/ 


(3.27) 


dt  \dqkJ  \dqkJ 

For  the  generalized  coordinates,  Sx ,  Sy ,  and  0,  the  Lagrangian  equations  of  motion  are  as 
follows: 

.a  m  o  t,  F\n  finr 

(3.28) 


X  o  X  3  m3  lo2  8y  8x 

m3  dx  -  2  ra3  u  dy - =  Q$x 
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m3  3y  +  2  m3  u)  5x  —  3  m3  lo2  5y  — 


9m^io25y2  3  m3  u)2  5x2 


2r0 


2  ra 


Qsy 


(3.29) 


mi 


6  +  3lu2  b 2  cos(0)  sin(0)  + 


3u2  fc3 
2  r0 


+  m2  (a2  9  +  3cu2  a2  cos(0)  sin(0)  —  3 % ^  °3  sin(0)^  =  Qe 


(3.30) 


The  generalized  forces,  Qfc;  are  found  according  to  Equation  3.31  in  order  to  complete 
the  equations  of  motion. 

Q*  =  Z$  ■§£  =  £*-  W  (3'31) 

i=i  i=i 

The  first  force  used  in  calculating  the  generalized  forces  is  the  force  of  drag  acting  on  each 
mass,  i. 


Fdragi  ~  0  C]ji  A{  p 


Vreh 


Vrelj 


(3.32) 


A  coefficient  of  drag,  Cd,  of  2  is  generally  the  value  used  in  analyzing  objects  in  space. 
Also,  the  density,  p,  is  approximated  using  an  atmospheric  model  written  by  David  Vallado 
at  the  Air  Force  Academy,  which  is  provided  in  Appendix  F.  The  relative  velocity,  Vrei,  is 
found  by  rotating  the  local  velocity  vector  to  the  inertial  reference  frame  and  subtracting 
the  velocity  of  the  Earth’s  atmosphere. 


Vreli 


Rot 


(u  t)  8ri  -  ( u Earth  z  x  Rot  {uj  t )  5r^j 


(3.33) 


The  rotation  matrix  between  the  inertial  and  body  frame  is: 


Rot  ((f)) 


f  —  sin(</>)  cos ((f))  0  \ 

cos  ((f))  sin(0)  0 

0  0  -1/ 


(3.34) 


The  drag  force  is  rotated  back  into  the  local  body  frame  by  multiplying  it  by  the  same 
rotation  matrix  shown  in  Equation  3.34. 

The  next  forces  considered  in  this  case  are  the  tension  forces  of  the  tether,  modelled 
as  spring  forces,  acting  on  each  mass.  These  forces  are  modelled  as  a  simple  spring  when 
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the  tether  is  taut  and  are  defined  as: 


spring  i 


Ti  El  (71  -  1) ,  71  >  1 

0,7i  <  1 


(3.35) 


spring 2 


f2  E2  (72  -  1) ,  72  >  1 
0, 72  <  1 

The  unit  vectors,  r,  along  which  the  forces  act  are: 


n 


1 

Ai 


b  sin (0)  —  Sx 
y  b  cos (9)  —  Sy  J 


(3.36) 


(3.37) 


1 


/ 

V 


—a  sin (0)  —  8x 
—a  cos (0)  —  8y 


(3.38) 


As  explained  in  Section  3.1,  the  moduli  of  elasticity  have  units  of  Newtons  per  millimeter 
squared,  but  the  cross-sectional  area  is  assumed  to  be  one  millimeter  squared.  Based  on 
this  assumption,  the  moduli  of  elasticity  are  in  units  of  Newtons  in  the  equations  of  motion 
and  in  the  computer  simulation  explained  in  Chapter  IV.  The  quantities,  7,  are  defined  by 
the  ratio  of  the  tether  length  at  a  certain  time  to  the  initial  tether  length: 


7i  = 


Ai_ 

A10 


72 


A2 

A20 


The  tether  lengths,  A,  at  a  certain  time  are  defined  by: 


(3.39) 

(3.40) 


Ai  =  yj {b  sin(0)  —  5x )2  +  (6  cos($)  —  Sy )2  (3-41) 

A2  =  \j (—a  sin(0)  —  8x )2  +  (—a  cos(0)  —  Sy )2  (3.42) 

These  forces  are  developed  from  equations  of  motion  derived  in  Dynamics  of  Space  Tethers 
by  Beletsky  and  Levin  for  the  motion  of  masses  at  the  ends  of  a  tether  (3:42-48).  These 
equations  are  altered  in  order  to  be  used  as  generalized  forces  using  the  chosen  generalized 
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coordinates.  The  motion  of  the  end  bodies  is  simplified  by  assuming  a  massless  tether.  This 
focuses  the  equations  of  motion  for  the  system  solely  on  the  end  masses.  Equations  3.35 
and  3.36  are  conditionals  which  make  an  allowance  for  a  slack  tether  since  compression  is 
not  a  factor  when  looking  at  the  tether  as  a  whole. 


Similarly,  a  dissipative  force  must  be  added  to  account  for  energy  losses  due  to  each 
tether.  This  energy  dissipation  allows  the  system  to  stabilize  over  time.  This  damping  is 
modelled  as: 


Fd 


amperi 


Cl  (f3  -  rij  ,71  >  1 
0, 7i  <  1 


(3.43) 


Fa 


,amper2 


C2  (r3  -  f2)  ,  72  >  1 

0, 72  <  1 


(3.44) 


with  damping  coefficients  for  tether  one  and  tether  two,  c\  and  C2,  respectively. 


Considering  the  forces  of  drag  and  tether  tensions  described  in  Equations  3.32  through 
3.44  and  substituting  into  Equation  3.31,  the  generalized  forces  are: 


Q  c  —  p  .  . _ — 

V ox  spnng\  „  • 

o  OX 


d  5rz 


j=>  u  ui  3  -tj 

T  damper 2  '  0  V  1“  drag 2 

a  ox 


(3.45) 


QSy  ~  Ft 


spring  i 


d  Sr  3 
d  5y 


+  F, 


spring 2 


9  V3  -  d8r3 

T  -t  damper  1 


+  Fdl 


ampere 


dSy 

dSr3 
d  5y 


+  Fd, 


■rag  3 


dSy 

d  8t3 
dSy 


(3.46) 


Qe  =  ~F( 


spnngi 


d8r\ 

d~6 


-F. 


spring 2 


dSr2  j=>  d5r  i 

T  -F  damper  1 


+  Fd, 


amper2 


d8r2 

~d¥ 


+  Fa 


■rag  1 


06 

0  8r\ 

~o¥ 


+  Fd 


rag  2 


06 
08r  2 

~of 


(3.47) 


The  above  generalized  forces  are  then  substituted  into  Equations  3.28  through  3.28  to  get 
the  final  equations  of  motion  for  the  two-dimensional  case  of  the  satellite  concept  depicted 
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in  Figure  3.1. 


m3  Sx  -  2  m3LuSy  -  3m3u;  &ySx  _  _I  <7^  p  Rot  t)  l/re;3fa  (3.48) 

r0  ^ 


((6  sin((9))  -  &c)  £1  (71  -  1)  Ax  *,71  >  1 


+ 


+ 


+ 


+ 


0, 71  <  1 

((-a  sin(0))  -  Sx)  E2  (72  -  1)  Ag  *,72  >  1 
0, 72  <  1 

ci  (Sx  —  io  Sy  —  bO  cos(0)  +  cob  cos  (0)  j  ,  71  >  1 
0, 71  <  1 

C2  (Sx  —  lu Sy  +  ad  cos(0)  —  lu  a  cos(0)^  ,  72  >  1 
0, 72  <  1 


••  •  2  9  m3  lu2  Sy2  3m3w2fa2 

m3Sy  +  Zm3uu  Sx  —  6m3u)  Sy - - - = 

2  rQ  2  rQ 

~  7  Cd3  Mp  WreA  Rot  (Wt)  Vrel36y 


+  < 


+ 


+ 


+ 


(( b  cos (0))  -  Sy)  Ei  (71  -  1)  X1 1,7i  >  1 
0, 71  <  1 

((-a  cos(d))  -  Sy)  E2  (72  -  1)  X21,72  >  1 
0,72  <  1 

ci  (Sy  +  lu  Sx  +  bO  sin(0)  —  lu  b  sin(0)  j  ,  71  >  1 
0, 71  <  1 

C2  (Sy  +  lu  Sx  —  a0  sin(0)  +  lu  a  sin(0)^  ,  72  >  1 
0, 72  <  1 


(3.49) 


/  O  Jl  7, 3  \ 

mi  [  62  0  +  3lu2  b 2  cos(0)  sin(0)  H - sin(0)  ) 

V  2r°  J 

+  m2  ( a2  6  +  2> lu2  a2  cos {6)  sin(0)  —  ^  sin^)^  = 


2r0 


(3.50) 


2^i 


Vrell 


Rot  (uut)  (vreli6x  b  cos(9)  -  Vreligy  b  sin(0)) 
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-  < 


-  \  CD2  A2  p  Vrel2  Rot  (u  t )  (Vrel2Sx  ~  a  COS (6)  +  Vrel2Sy  a  Sm(0)) 

[(6  sin(0)  —  Sx)  b  cos (0)  —  ( b  cos (0)  —  Sy)  b  sin(0)]  E\  (71  —  1)  A^qi  >  1 

0, 71  <  1 

[(—a  sin(0)  —  Sx)  —  a  cos(0)  +  (—a  cos(0)  —  Sy)  a  sin(0)]  E2  (72  —  1)  A^"1,72  >  1 

0, 72  <  1 

ci  (b  cos (0)  (fix  —  uj  Sy  —  b9  cos(0)  +  uj b  cos (fl)'j 'j  , 71  >  1 

0,7!  <1 

ci  (^—b  sin(0)  [Sy  +  uj  Sx  +  b  6  sin(0)  —  uj  b  sin(0)^  j  ,  71  >  1 

0,7i  <  1 

C2  a  cos (0)  (Sx  —  uj  Sy  —  b9  cos(0)  +  uj  b  cos (fl)'j 'j  ,  72  >  1 

0, 72  <  1 

C2  (a  sin(0)  (Sy  +  uj  Sx  +  b0  sin(0)  —  ujb  sin(0)^  ,  72  >  1 


+  < 


+  < 


+ 


+ 


0, 72  <  1 


These  equations  of  motion  are  used  in  a  fourth-order  Runge-Kutta  numerical  inte¬ 
grator  which  propagates  the  system  parameters,  Sx,  Sy,  and  6  forward  in  time.  These 
equations  of  motion  are  solved  for  the  accelerations  before  being  used  in  the  numerical 
integrator,  which  is  described  further  in  Section  4.3. 
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IV.  Methodology 

The  program  used  to  simulate  the  attitude  of  this  satellite  concept  is  described  in  this  sec¬ 
tion.  While  defining  the  equations  of  motion  is  paramount  to  creating  a  simulation,  writing 
this  program  has  been  the  pith  of  this  study.  Contained  in  this  chapter  is  a  description 
of  the  simulation  program  in  terms  of  the  satellite  values  tested.  A  general  program  algo¬ 
rithm  is  provided  that  describes  the  simulation  at  the  conceptual  level.  The  chapter  also 
describes  the  approach  taken  to  run  the  simulations  whose  results  are  found  in  Chapter  V 
and  discusses  the  method  of  program  validation.  The  actual  code  of  the  program  is  provided 
in  Appendices  A  through  H.  The  program  code  has  been  separated  into  each  subprogram 
and  is  generally  listed  in  the  order  in  which  it  is  called.  Contained  within  each  program 
code  in  the  appendices  is  an  algorithm  as  well  as  a  list  of  variables,  constants,  and  coupled 
programs. 

4.1  Satellite  Characteristics 

The  characteristics  of  the  satellites  used  in  this  simulation  represent  modest  estima¬ 
tions  of  microsatellites  with  a  gravity  gradient  boom.  Since  no  satellite  similar  to  that 
depicted  in  Figure  3.1  exists,  each  variable  either  was  defined  by  comparison  to  similar 
satellite  components  or  by  determining  a  range  of  values  through  iterative  simulation.  The 
satellite  specifications  used  in  this  study  are  listed  in  Table  4.1. 

The  listed  units  are  the  same  units  used  within  the  simulation  program  with  the 
exception  of  the  moduli  of  elasticity.  As  stated  in  the  assumptions,  the  listed  moduli  of 
elasticity  are  first  multiplied  by  the  assumed  cross-sectional  area  of  the  tethers,  which  is  one 
millimeter  squared,  before  being  used  in  the  simulation.  Also,  the  listed  moduli  of  elasticity 
and  damping  coefficients  are  nominal  values  for  each  satellite.  These  were  found  by  iteration 
of  the  simulation,  which  is  described  further  in  Chapter  V.  For  these  characteristics,  many 
values  are  possible  within  a  range  dependent  on  stability  and  performance  specifications. 
This  is  also  discussed  in  Chapter  V. 

The  first  satellite,  Satellite  I  or  Satl,  to  be  modelled  represents  a  very  general  case. 
The  values  for  this  satellite  are  a  compromise  between  reality  and  simplicity.  For  example, 
these  general  values  make  the  calculations  for  the  position  of  the  center  of  mass  (COM)  and 
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Table  4.1  Satellite  Specifications 


Variable 

Symbol 

Satellite  I 

Satellite  II 

DumbSat 

Units 

Tip-Mass 

mone 

10 

3 

50 

kg 

Main  Bus  Mass 

mtwo 

100 

70 

50 

kg 

Balloon  Mass 

mthree 

1 

0.5 

1 

kg 

Tip-Mass  Frontal  Area 

areaone 

0.1 

0.0231 

0.2 

9 

m 

Main  Bus  Frontal  Area 

areatwo 

0.5 

0.2275 

0.2 

m2 

Balloon  Frontal  Area 

areathree 

5 

2 

2 

9 

m 

Boom  Length 

a  +  b 

6 

5 

5 

m 

Tip  to  COM  Distance 

b 

5.5 

4.7945 

2.5 

m 

Bus  to  COM  Distance 

a 

0.5 

0.2055 

2.5 

m 

T1  Taut  Length 

lonenote 

6 

5 

5 

m 

T2  Taut  Length 

ltwonot 

6 

5 

5 

m 

T1  Modulus  of  Elasticity 

Eone 

0.002065 

0.0003 

0.00317 

N /mm2 

T2  Modulus  of  Elasticity 

Etwo 

0.0202 

0.006 

0.00317 

N /mm2 

T1  Damping  Coefficient 

cone 

0.0007 

0.0007 

0.0016 

Ns/m 

T2  Damping  Coefficient 

etwo 

0.0007 

0.0007 

0.0016 

Ns/m 

inertia  matrix  easier.  Although  the  values  are  realistic,  they  do  not  represent  a  microsatellite 
as  well  as  the  second  satellite.  The  main  criteria  for  the  first  satellite  is  that  it  be  a  stable 
gravity  gradient  satellite.  The  criteria  for  stability  is  to  have  the  principle  moments  of 
inertia  be  such  that  the  major  axis  (C)  points  normal  to  the  orbit,  the  intermediate  axis 
(B)  points  in  the  velocity  direction,  and  the  minor  axis  (A)  points  in  the  radial  direction. 
This  stability  criteria  is  explained  further  in  Section  2.2  and  Section  4.4.  The  equations  and 
values  for  the  moments  of  inertia  are  outlined  in  Section  4.4. 

Satellite  II,  or  Satll,  uses  more  realistic  values  to  increase  the  validity  of  the  results. 
For  example,  the  mass  of  the  main  satellite  bus  was  taken  from  a  range  of  values  given  by 
Surrey  Satellite  Technology  Ltd.  Their  information  was  used  because  their  flight-proven 
systems  are  designed  with  the  microsatellite  concept  in  mind  (4).  In  addition,  the  speci¬ 
fications  for  the  gravity  gradient  boom  were  defined  based  on  information  provided  on  an 
actual  boom  created  by  Northrop- Grumman  specifically  for  microsatellites  (11). 

The  last  satellite  tested,  DumbSat,  is  essentially  a  dumbbell  configuration.  This  satel¬ 
lite  tests  a  completely  different  configuration  from  Satellite  I,  adding  a  breadth  to  the 
simulation  results.  Also,  the  tether  characteristics  for  a  stable  configuration  are  more  intu¬ 
itive.  Since  both  ends  of  DumbSat  are  equal  in  mass,  length,  and  frontal  area,  the  attached 
equal  length  tethers  should  have  equal  characteristics. 
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4-2  Simulation  Approach 

In  the  simulations  run,  each  satellite  described  in  Table  4.1  is  simulated  with  the  initial 
generalized  coordinates  in  Table  4.2.  The  velocity  terms  for  each  generalized  coordinate  are 


Table  4.2  Initial  Generalized  Coordinates 


Generalized  Coordinate 

Variable 

Satl 

Satll 

Dumb  Sat 

Units 

Balloon  Position  in  x-direction 

5x 

-3.25 

-2.6028 

-3.75 

meters 

Balloon  Position  in  y-direction 

Sy 

4.7631 

4.1522 

2.1651 

meters 

Pitch  of  Gravity  Gradient  Satellite 

e 

7t/6 

7 r/6 

7r/6 

radians 

assumed  to  be  zero  in  all  cases.  Each  simulated  satellite  also  begins  at  the  same  position 
as  described  in  Table  4.3  for  the  satellite  center  of  mass.  These  initial  coordinates  simulate 


Table  4.3  Initial  Inertial 

Coordinates 

Inertial  Coordinate 

Value 

Units 

Rx 

6578.1 

kilometers 

Ry 

0 

kilometers 

Rz 

0 

kilometers 

vx 

0 

kilometers/sec 

Vy 

7.7843 

kilometers/sec 

vz 

0 

kilometers/sec 

a  circular,  equatorial  satellite.  Although  an  equatorial  orbit  is  not  an  assumption,  a  non- 
equatorial  orbit  would  provide  the  same  results  since  a  spherical  Earth  is  assumed.  In 
addition,  the  number  of  steps  and  step  size  varied  depending  on  what  was  the  aim  of  the 
simulation.  For  example,  to  get  a  general  picture  of  a  certain  case,  fewer  steps  allows  the 
simulation  to  run  faster.  On  the  other  hand,  more  steps  are  required  when  measurements 
are  taken  directly  from  the  graph.  The  number  of  steps  is  also  constrained  by  the  numerical 
integrator,  discussed  further  in  Section  4.3.  A  small  number  of  steps  decreases  accuracy 
because  there  are  not  enough  data  points  to  produce  an  accurate  graph,  but  a  large  number 
of  steps  results  in  a  longer  run  time  for  the  simulation.  The  number  of  steps  and  interval 
for  each  presented  simulation  result  is  detailed  in  the  corresponding  section  of  the  report  or 
graph. 
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4.3  Program  Algorithm 

The  algorithm  of  this  program  stems  from  a  basic  orbit  prediction  program  self- written 
for  an  undergraduate  degree.  The  original  program  uses  a  method  called  Cowell  to  predict 
and  track  the  position  of  a  satellite  using  the  two-body  equations  of  motion  and  certain 
orbital  perturbations,  like  drag  and  the  effects  of  the  Earth’s  oblateness.  That  algorithm 
has  been  adapted  to  use  the  equations  of  motion  described  in  Section  3.4  and  produce  results 
that  are  used  to  determine  this  satellite  concept’s  feasibility. 

The  numerical  integrator  is  to  simulate  the  motion  of  the  satellite  was  coded  in  Matlab. 
The  equations  of  motion  derived  in  Section  3.4  are  rearranged  slightly  to  solve  for  acceler¬ 
ation  and  placed  into  a  fourth  order  Runge-Kutta  numerical  integrator.  The  Runge-Kutta 
numerical  integrator  is  based  on  Equation  4.1  (2:414). 

Xn+i  =  Xn  +  —  ( k\  +  2  +  2  k%  +  Aq)  (4-1) 

In  this  equation,  the  k’s  are  defined  as: 


k1  =  hf(tn,Xn)  (4.2) 

^2  —  h  f  ^tn  +  — ,  Xn  +  — ^  (4-3) 

fc3  =  /*/(*n  +  ^,Xn  +  ^) 

Aq  =  h  f  (tn  +  /i,  Xn  +  ks) 


where  h  is  the  step  size  and  the  function  /  (£,  X)  is  equal  to  the  derivative  of  the  state 
vector  or  X.  This  integrator  is  normally  used  for  a  state  vector,  X,  with  six  variables  - 
three  position  vectors  and  three  rate  vectors.  According  to  Bate,  Mueller,  and  White,  the 
Runge-Kutta  algorithm  is  stable  and  has  a  small  truncation  error  (2:415). 

To  use  this  integrator,  the  simulation  program  begins  by  declaring  all  constants  and 
satellite  variables  as  listed  in  Table  4.1.  Second,  the  number  of  steps  and  the  step  size  are 
defined  based  on  the  period  of  the  orbit  and  the  number  of  orbits  intended  for  inclusion. 
Then,  the  program  iterates  through  the  number  of  steps  in  increments  of  the  defined  step 
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size,  calling  subprograms  that  determine  the  future  position  of  the  center  of  mass  as  well  as 
the  three  masses  of  the  satellite  system.  These  subprograms  calculate  the  state  vector  and 
its  derivative,  which  depends  on  the  position  of  the  balloon.  If  the  balloon’s  position  causes 
either  tether  to  be  taut,  then  that  tether’s  spring  and  damping  forces  are  included  in  the 
equations  of  motion.  The  state  vector  and  its  derivative  are  then  used  in  Equation  4.1  to 
calculate  the  future  state  vector.  At  each  time  step,  the  state  vector  is  recorded  so  that  at 
the  end  of  the  program  a  plot  can  be  made  of  the  progress  of  Sx,  5y ,  and  9. 


4-4  Program  Validation 

In  order  to  determine  whether  the  Matlab  simulation  program  functions  properly,  a 
case  has  been  run  which  has  known  results.  The  case  of  a  satellite  with  a  gravity  gradient 
boom  without  the  effects  of  drag  has  a  known  response  to  initial  conditions  so  it  was 
used  to  test  the  computer  simulation  of  this  satellite  concept.  For  this  case,  the  forces 
due  to  the  tether  have  been  neglected  regardless  of  the  position  of  the  balloon  mass.  The 
resulting  oscillations  in  pitch  angle,  0,  have  a  predictable  frequency  that  is  dependent  on 
the  principal  moments  of  inertia.  Using  the  specifications  for  each  satellite  listed  in  Table 
4.1  and  assuming  a  box-shaped  main  satellite  bus  and  tip-mass,  the  principal  moments  of 
inertia  were  calculated  using  the  following  equation: 


Ixx  —  lyy  —  Izz  —  22  m 


(4.4) 


where  Area  is  the  frontal  area  (1:702).  The  resulting  principal  moments  of  inertia  for  the 
gravity  gradient  portions  of  each  satellite  are  listed  in  Table  4.4.  These  values  are  based  on 
6  equal  to  zero.  The  frequency  of  oscillations  in  the  pitch  direction,  np,  has  a  predictable 


Table  4.4  Principal  Moments  of  Inertia 


Principal  Axis 

Satellite  I 

Satellite  II 

DumbSat 

Units 

Minor  Axis  (B) 

336 

74.583 

628 

kg/m2 

Intermediate  Axis  (A) 

8.5 

2.666 

3.333 

kg/m2 

Major  Axis  (C) 

336 

74.583 

628 

kg/m2 
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value  as  defined  by  Equation  4.5  and  has  units  of  radians  per  second  (17:153). 


nn 
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The  resulting  predicted  oscillation  frequencies  for  the  pitch  angle  are  listed  in  Table  4.5. 


Table  4.5  Frequencies  of  Pitch  Angle  Oscillation 


Result 

Satellite  I 

Satellite  II 

Dumb  Sat 

Units 

Predicted 
Simulated 
Percent  Error 

0.002024 

0.00190978 

5.68 

0.002013 

0.00190851 

5.17 

0.002044 

0.00190974 

6.60 

rad/sec 

rad/sec 

Percent 

The  simulated  values  were  calculated  by  examining  the  output  graph  of  pitch  angle, 
0,  versus  time.  The  output  for  Satellite  I  is  provided  in  Figure  4.1,  and  the  simulated 
oscillation  frequency  is  listed  in  Table  4.5.  The  interval  in  this  test  is  two  periods  with  400 
steps.  Sequential  maximums  are  recorded  to  determine  the  period  in  terms  of  number  of 


0  vs.  Time  Step  for  Satl 


Figure  4.1  Simulated  Results  for  Oscillation  Frequency  Validation  for  Satl 
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steps,  X.  The  oscillation  frequency  is  then  found  by  Equation  4.6. 


nv  ( simulated )  =  — — ^ — —  (4-6) 

The  step  size  has  units  of  seconds  per  step,  and  the  resulting  frequency  from  Equation  4.6 
is  in  units  of  radians  per  second. 

This  test  case  was  also  run  for  Satellite  II.  The  results  are  also  in  Tables  4.4  and  4.5, 
while  the  output  for  Satellite  II  is  provided  in  Figure  4.2.  This  test  case  was  also  run  for 


0  vs.  Time  Step  for  Satll 


Figure  4.2  Simulated  Results  for  Oscillation  Frequency  Validation  for  Satll 
DumbSat.  The  results  are  also  in  Tables  4.4  and  4.5.  The  output  is  provided  in  Figure  4.3. 

It  is  important  to  note  that  all  of  these  graphs  show  pitch  angle  oscillations  varying 
between  tv/ 6  and  negative  7r/6,  and  these  oscillations  do  not  dampen.  Hughes,  mentioned 
previously  in  Section  2.2,  points  out  that  certain  satellite  configurations  produce  oscillations 
with  marginal  stability.  One  of  these  configurations  is  when  k 3  is  approximately  zero  and 
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0  vs.  Time  Step  for  DumbSat 


Figure  4.3  Simulated  Results  for  Oscillation  Frequency  Validation  for  DumbSat 

when  k\  is  positive  (6:302).  This  is  equivalent  to  a  slender  rod  with  the  long  axis  pointed  in 
the  radial  direction.  As  seen  in  Figure  3.1,  these  satellites  closely  resemble  slender  rods  with 
the  long  axis  pointed  in  the  radial  direction.  Furthermore,  marginal  stability  is  characterized 
by  an  inability  to  dampen  oscillations  in  a  specified  direction.  As  shown  in  the  results  of 
this  test  case,  none  of  the  satellites  have  pitch  angle  graphs  that  decrease  in  amplitude, 
demonstrating  marginal  stability. 

These  three  tests  demonstrate  the  accuracy  of  the  system  in  that  they  have  very  low 
percent  error.  Using  the  three  satellites  also  shows  consistency  and  a  breadth  in  usability 
for  the  simulation.  The  error  in  the  simulated  results  stems  from  the  assumptions  made, 
mainly  that  there  is  no  motion  outside  of  the  orbital  plane.  Because  of  this  assumption, 
the  dynamics  and  simulation  do  not  account  for  any  rolling  or  yawing  motions.  Also,  these 
satellites  are  quite  similar  and  therefore  have  similar  oscillation  frequencies. 
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V.  Results 


After  the  program  was  written  and  validated,  the  focus  was  placed  on  the  main  research 
objective,  determining  whether  passive  attitude  stabilization  can  be  achieved  by  manipu¬ 
lating  aerodynamic  and  gravity  gradient  torques.  Attention  was  first  placed  on  determining 
the  effect  of  altering  the  modulus  of  elasticity  for  the  two  tethers  in  the  satellite  system. 
Once  a  stable  oscillation,  or  acceptable  steady  state  pitch  angle,  was  observed,  the  damping 
coefficients  were  then  altered  to  determine  their  impact  on  the  system.  It  was  also  im¬ 
portant  to  ascertain  a  legitimate  damping  of  the  pitch  angle  over  a  limited  time,  which  is 
the  same  as  an  acceptable  settling  time.  These  factors  were  then  broken  down  into  forces 
applied  on  each  mass  for  each  generalized  coordinate.  Once  the  influence  of  these  forces 
was  understood,  the  regions  of  stability  for  the  moduli  of  elasticity  and  damping  coefficients 
were  determined  for  the  three  satellites.  This  analysis  was  first  performed  on  Satellite  I 
so  the  first  set  of  results  discussed  are  all  for  this  satellite.  The  results  of  Satellite  II  and 
DumbSat  are  provided  following  the  analysis  of  Satellite  I. 

5. 1  Modulus  of  Elasticity 

The  modulus  of  elasticity  was  the  first  variable  to  be  inspected.  This  is  the  first 
simulation  to  include  the  effects  of  the  attached  balloon  mass  and  air  drag  acting  on  all 
three  masses.  The  damping  coefficient  was  held  constant  at  zero  to  isolate  the  effects  of  the 
moduli  of  elasticity  on  the  system.  This  simulation  was  iterated  until  an  acceptable  value 
for  the  moduli  of  elasticity  was  found  which  produced  a  stable  system.  For  these  iterative 
simulations,  the  inspected  time  interval  was  five  periods  with  1500  steps  and  a  step  size 
of  17.7  seconds.  After  iterating  through  this  case,  an  acceptable  pitch  angle  output  was 
established  for  E\  equal  to  0.002065  N /mm2  and  E2  equal  to  0.0202  N /mm2,  shown  in  the 
output  graph  of  pitch  angle, 0,  in  Figure  5.1. 

The  nominal  moduli  of  elasticity  found  are  very  small  values  when  compared  to  typical 
materials  used  for  space  tethers.  As  shown  in  Table  5.1,  typical  moduli  of  elasticity  have 
values  several  orders  of  magnitude  greater  than  the  values  found  for  the  nominal  case.  This 
is  due  in  part  on  the  application  of  typical  space  tethers.  In  general,  space  tethers  are  not 
intended  to  be  used  as  springs  but  are  used  more  to  hold  components  together  at  a  certain 
distance.  In  this  satellite  concept,  the  moduli  of  elasticity  are  intended  to  act  as  springs 
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6  vs.  Time  Step  for  Satl 


Time  Steps  (x17.7sec/step) 

Figure  5.1  Nominal  Moduli  of  Elasticity  for  Satl 


Table  5.1  Moduli  of  Elasticity  for  Space  Tethers  (3:35) 


Material 

Modulus  of  Elasticity  (kN /mm2) 

Percent  Elongation 

Kevlar-49 

130 

2.5 

Alloyed  Aluminum 

70 

10 

Stainless  Steel 

200 

1.4 

so  that  a  damping  coefficient  can  be  added  in  a  realistic  manner.  Table  5.2  lists  moduli 
of  elasticity  closer  to  those  found  in  the  simulation.  Although  there  is  no  known  flight 


Table  5.2  Moduli  of  Elasticity  for  Possible  Tethers 

(10) 

Material 

Modulus  of  Elasticity  (N /mm2) 

Percent  Elongation 

Silastic  (R)  24005  Silicon  Rubber 

0.0588 

950 

Silastic  (R)  24064  Silicon  Rubber 

0.1449 

790 

Silastic  (R)  29051  Silicon  Rubber 

0.2373 

436 

experience  of  these  materials,  they  are  able  to  operate  in  temperatures  well  below  freezing 

(10). 

The  average  of  the  pitch  angles  for  the  entire  interval  is  0.0412  radians  or  2.36  degrees 
and  is  within  a  very  small  margin  of  error  from  the  ideal  zero  degree  angle.  This  average 
pitch  angle  correlates  to  a  gravity  gradient  satellite  with  the  main  satellite  in  a  nadir  pointing 
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configuration.  This  system  is  acceptable  since  it  is  marginally  stable  as  the  oscillations  do 
not  extend  beyond  the  initial  amplitude.  Also  in  Figure  5.1,  the  pitch  angle  is  oscillating 
between  positive  7r/6,  which  is  the  initial  pitch  angle  value,  and  negative  7r/6.  This  correlates 
to  the  gravity  gradient  portion  of  the  satellite  oscillating  like  an  undamped  pendulum.  This 
is  expected  for  a  marginally  stable  gravity  gradient  satellite. 

Figure  5.2  shows  the  motion  of  the  balloon  mass  for  the  coordinate  Sx.  The  values  of 


5  x  vs.  Time  Steps 


Figure  5.2  Nominal  Moduli  of  Elasticity  for  Satl 

Sx  continue  to  oscillate  in  this  case  because  the  pitch  angle  of  the  gravity  gradient  portion 
of  the  satellite  oscillates  in  an  undamped  manner.  Figure  5.2  shows  that,  while  the  taut 
length  of  the  tether  is  set  to  six  meters,  the  steady  state  length  is  about  25.25  meters.  This 
is  as  if  the  spring  has  reached  its  maximum  length  or  is  completely  stretched.  Table  5.1 
also  includes  percent  elongation  for  typical  tether  materials  while  Table  5.2  includes  values 
for  possible  materials.  The  main  difference  in  values  is  again  attributed  to  the  difference  in 
applications.  The  tethers  in  this  satellite  concept  are  intended  to  act  as  springs,  and  springs 
typically  extend  much  further  than  the  initial  taut  length. 
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Figure  5.3  shows  the  motion  of  the  balloon  for  Sy  over  the  simulated  interval.  The 
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Figure  5.3  Nominal  Moduli  of  Elasticity  for  Satl 

balloon  mass  is  oscillating  in  the  vertical  direction  between  approximately  positive  and 
negative  two  meters.  This  is  also  due  to  the  oscillation  in  the  pitch  angle. 

An  analysis  was  also  conducted  for  unacceptable  values  for  moduli  of  elasticity.  When 
the  modulus  of  elasticity  was  well  out  of  range,  it  was  easy  to  see  that  the  pitch  angle  spun 
off  into  infinity.  An  example  is  shown  in  Figure  5.4  where  E\  is  equal  to  0.0030975  N /mm2 
and  E2  is  equal  to  0.0202  N /mm2.  The  pitch  angle  values  are  limited  to  between  negative 
two  pi  and  positive  two  pi  in  the  figure  so  when  the  satellite  tumbles  past  two  pi  radians 
or  360  degrees,  the  graph  shows  pitch  angles  going  from  two  pi  directly  back  to  zero.  Since 
Figure  5.4  shows  the  pitch  angle  continually  doing  this,  the  simulated  satellite  is  actually 
unstable  and  is  tumbling  end-over-end.  This  shows  the  volatility  of  the  system  since  the 
modulus  of  elasticity  for  tether  one  is  only  fifty  percent  more  than  its  nominal  value  of 
0.002065  N /mm2. 
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Figure  5.4  Modulus  of  Elasticity  Out  of  Range  by  50  Percent  for  Satl 

Furthermore,  it  was  necessary  to  observe  the  average  of  the  pitch  angles  to  ensure  that 
it  was  stabilizing  to  zero  degree.  This  corresponds  to  the  nadir  pointing  orientation.  An 
example  of  pitch  angles  oscillating  about  a  non-zero  value  is  shown  in  Figure  5.5  where  E\ 
is  equal  to  0.0022715  N /mm2  and  E2  is  equal  to  0.0202  N /mm2.  The  average  in  this  case  is 
-1.5504  radians,  which  corresponds  to  oscillations  about  ninety  degrees  off  the  vertical  axis. 
If  damping  were  added  to  this  case,  the  satellite  would  be  oriented  with  the  main  satellite 
mass  pointed  in  the  ram  direction  and  the  tip  mass  trailing  directly  behind  in  the  horizontal 
direction.  This  may  be  desired  in  certain  missions,  but  that  is  not  the  intended  steady  state 
orientation  of  this  simulation. 

The  impact  of  increasing  and  decreasing  either  moduli  of  elasticity  is  also  important 
to  know  when  attempting  to  alter  the  system  to  meet  certain  criteria.  As  shown  in  Figure 
5.5,  increasing  the  first  tether’s  modulus  of  elasticity  decreases  the  average  pitch  angle; 
therefore,  a  decrease  in  the  first  tether’s  modulus  of  elasticity  increases  the  average  pitch 
angle.  The  opposite  is  true  for  the  second  tether’s  modulus  of  elasticity.  A  change  in  the 
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Figure  5.5  Modulus  of  Elasticity  Out  of  Range  by  10  Percent  for  Satl 

second  modulus  of  elasticity  has  a  directly  proportional  impact  on  the  average  pitch  angle. 
For  example,  increasing  the  second  modulus  of  elasticity  increases  the  average  pitch  angle. 

5.2  Damping  Coefficient 

The  next  variables  introduced  into  the  system  are  the  damping  coefficients  for  both 
tethers,  which  simulate  energy  dissipation  in  the  tethers  due  to  real-world  imperfections. 
This  first-level  inspection  of  damping  coefficient  is  important  to  establish  that  this  satellite 
characteristic  is  the  major  factor  in  damping  down  the  oscillating  pitch  angle.  The  criteria 
for  this  variable  is  that  the  average  pitch  angle  stays  close  to  zero  degree.  Second,  the 
settling  time  for  the  satellite  system  must  be  within  three  days,  which  equates  to  about  fifty 
orbits  for  a  200  kilometer  altitude,  circular  orbit.  The  last  criteria  is  that  the  forces  exerted 
on  each  mass  and  the  value  for  the  damping  coefficient  cannot  be  too  excessive.  In  terms  of 
the  former,  large  forces  exerted  on  the  masses  could  cause  the  tethers,  tether  connections, 
or  on-board  instruments  to  fail.  In  terms  of  the  latter,  the  damping  coefficient  should  not 
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be  too  large  as  this  variable  is  only  meant  to  model  imperfections  in  the  tether  material. 
In  reality,  energy  would  be  lost  through  the  tether  in  the  form  of  radiation  to  space  so  an 
overly  large  damping  coefficient  would  be  unrealistic. 

During  these  simulations,  the  inspected  time  interval  was  fifty  periods  with  10,000 
steps  and  a  step  size  of  26.55  seconds.  This  case  was  run  with  the  nominal  moduli  of 
elasticity  shown  in  Table  4.1  and  depicted  in  Figures  5.1  through  5.3.  Nominal  values  for 
damping  coefficients  were  found  around  0.0007  Ns/m  for  both  tethers  for  Satellite  I.  The 
results,  shown  in  Figure  5.6,  had  an  average  of  0.02  radians  and  settled  to  0.0219  radians 
after  a  little  over  two  days. 
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Figure  5.6  Stable  Damping  Coefficients  for  Satl 


Forces  on  m  sin  X-Direction 


Figure  5.7  Forces  Acting  on  Balloon  Mass  in  X-Direction  for  Satl 


The  forces  acting  on  each  generalized  coordinate  are  shown  in  Figures  5.7  through 
5.12.  The  first  graph  shows  the  impacts  of  each  force  on  6x.  The  maximum  force  felt  by  the 
balloon  in  the  x-direction  is  about  two  Newtons.  The  major  forces  acting  on  the  balloon 
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Figure  5.8  Zoomed  In  View  of  Forces  Acting  on  Balloon  Mass  in  X-Direction  for  Satl 


mass  in  the  x-direction  are  drag  and  the  spring  force  of  tether  one.  The  drag  force  pulls 
the  balloon  mass  in  the  negative  velocity  direction  while  the  spring  force  from  tether  one 
pulls  it  in  the  positive  velocity  direction.  The  forces  due  to  the  second  spring,  both  tether 
damping  effects,  and  two-body  effects  are  very  small  so  that  region  is  expanded  in  Figure 
5.8. 

Figure  5.9  shows  the  impact  of  each  force  on  the  balloon  mass  in  the  Sy  coordinate. 
The  maximum  force  felt  by  the  balloon  in  the  y-direction  is  about  1.2  Newtons.  The  major 
forces  acting  on  the  balloon  mass  in  the  y-direction  are  the  damping  forces  of  each  tether. 
This  shows  that  the  damping  forces  of  each  tether  are  the  dominating  mechanism  in  settling 
the  balloon  mass  to  a  certain  value  in  the  radial  or  vertical  direction.  Each  force  except  the 
first  and  second  dampers  are  very  small  so  the  region  close  to  zero  Newtons  is  expanded 
in  Figure  5.10.  These  show  that  the  damping  force  and  spring  force  for  each  tether  act  in 
opposition,  as  is  the  case  for  spring-damper  systems. 

Figure  5.11  shows  the  impacts  of  each  force  on  6.  The  maximum  force  felt  by  the 
gravity  gradient  portion  of  the  satellite  system  is  102  Newtons.  The  dominating  forces  for 
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Figure  5.9  Forces  Acting  on  Balloon  Mass  in  Y-Direction  for  Satl 
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Figure  5.10  Zoomed  In  View  of  Forces  Acting  on  Balloon  Mass  in  Y-Direction  for  Satl 
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Forces  acting  on  m  land  m  Jn  Pitch  Angle  Direction 
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Figure  5.11  Forces  Acting  on  Tip-Mass  and  Main  Bus  Mass  along  6  for  Satl 

this  coordinate  are  the  spring  forces  of  both  tethers.  Each  force  except  the  first  and  second 
springs  are  very  small  so  that  region  is  expanded  in  Figure  5.12. 

In  analyzing  the  effect  of  damping  coefficient  on  pitch  angle  stability,  it  can  be  seen 
that  there  are  upper  and  lower  ranges.  The  lower  range  fails  to  meet  the  settling  time 
requirement  of  three  days.  This  lower  limit  is  dependent  on  other  satellite  characteristics, 
like  moduli  of  elasticity,  the  masses  of  each  system  part,  and  the  gravity  gradient  boom 
length.  In  the  case  of  Satellite  I  and  using  the  nominal  moduli  of  elasticity  in  Table  4.1,  the 
lower  limits  for  the  damping  coefficients  are  illustrated  in  Figure  5.13. 
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Figure  5.12  Zoomed  In  View  of  Forces  Acting  on  Tip-Mass  and  Main  Bus  Mass  along  6 
for  Satl 
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Figure  5.13  Stability  Region  for  Damping  Coefficients  for  Satl 
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On  the  other  hand,  the  upper  range  produces  pitch  angles  which  are  above  the  de¬ 
sired  average.  When  the  nominal  damping  coefficients  are  multiplied  by  ten,  the  result 
is  as  shown  in  Figure  5.14.  The  steady  state  pitch  angle  for  this  case  is  3.18  or  about  n 
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Figure  5.14  Damping  Coefficients  Out  of  the  Desired  Range  for  Satl 

radians.  This  is  the  same  as  the  satellite  being  oriented  upside-down  at  steady  state,  which 
is  possible  in  normal  gravity  gradient  satellites  as  well.  One  problem  with  normal  gravity 
gradient  satellites  is  that  the  upside-down  orientation  is  actually  an  equilibrium  point  and 
has  been  the  result  for  some  gravity  gradient  satellites.  If  the  damping  coefficient  is  not 
modelled  correctly,  an  upside-down  orientation  is  a  possibility.  Just  as  normal  gravity  gra¬ 
dient  satellites  require  a  mechanism  to  avoid  the  upside-down  orientation,  gravity  gradient 
satellites  with  a  tethered  balloon  may  require  such  a  mechanism  if  the  tether  characteristics 
are  not  modelled  and  designed  correctly. 

In  addition,  the  forces  applied  to  each  mass  were  also  increased.  The  maximum  force 
felt  by  the  balloon  in  the  x-direction  is  about  1.9  Newtons,  as  shown  in  Figure  5.15.  The 
dominating  forces  on  the  balloon  in  this  direction  is  the  spring  force  of  the  first  tether  and 
the  drag  force,  which  remains  the  same  from  the  nominal  damping  coefficients  results.  Each 
force  except  the  first  spring  force  and  drag  are  very  small  so  that  region  is  expanded  in 
Figure  5.16. 
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Figure  5.15  Forces  Acting  on  Balloon  Mass  in  X-Direction  For  Damping  Coefficient  Out 
of  Range  for  Satl 

Figure  5.17  shows  the  impacts  of  each  force  on  Sy.  The  maximum  force  felt  by  the 
balloon  in  the  y-direction  is  1.1  Newtons.  The  major  contributing  forces  in  this  direction  are 
the  damping  forces  from  each  tether,  as  is  the  case  from  the  nominal  damping  coefficients 
results.  Each  force  except  the  first  and  second  dampers  are  very  small  so  that  region  is 
expanded  in  Figure  5.18. 

Figure  5.19  shows  the  impacts  of  each  force  on  the  pitch  angle.  The  maximum  force 
felt  by  the  gravity  gradient  portion  of  the  satellite  system  is  157  Newtons.  Just  like  the 
results  for  the  nominal  damping  coefficients,  the  spring  forces  from  both  tethers  are  the 
dominating  forces.  Each  force,  except  the  first  and  second  springs,  are  very  small  so  that 
region  is  expanded  in  Figure  5.20.  As  shown,  the  satellite  masses  experience  a  maximum 
force  greater  by  a  factor  of  fifty  percent  when  the  damping  coefficients  were  increased  by  an 
order  of  magnitude. 
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Figure  5.16  Zoomed  in  View  of  Forces  Acting  on  Balloon  Mass  in  X-Direction  For  Damp¬ 
ing  Coefficient  Out  of  Range  for  Satl 
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Figure  5.17  Forces  Acting  on  Balloon  Mass  in  Y-Direction  For  Damping  Coefficient  Out 
of  Range  for  Satl 
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Forces  on  m  3in  Y-Direction 
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Figure  5.18  Zoomed  in  View  of  Forces  Acting  on  Balloon  Mass  in  Y-Direction  For  Damp¬ 
ing  Coefficient  Out  of  Range  for  Satl 
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Figure  5.19  Forces  Acting  on  Tip-Mass  and  Main  Bus  Mass  For  Damping  Coefficient  Out 
of  Range  for  Satl 
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Figure  5.20  Zoomed  in  View  of  Forces  Acting  on  Tip-Mass  and  Main  Bus  Mass  For  Damp¬ 
ing  Coefficient  Out  of  Range  for  Satl 
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5.3  Stability  Region 

Although  values  for  the  moduli  of  elasticity  and  damping  coefficients  are  provided  in 
Table  4.1,  these  are  not  the  only  values  which  work.  The  values  used  in  the  initial  program 
validation  were  chosen  because  they  lie  in  the  middle  of  the  stability  region.  In  fact,  a  range 
of  values  had  to  be  established  through  iterative  simulation.  By  running  the  simulation  for 
many  different  tether  characteristics,  a  region  of  stability  was  determined.  The  criteria  for 
this  stability  region  is  based  on  the  steady  state  pitch  angle  and  the  settling  time. 

The  criteria  established  for  the  modulus  of  elasticity  uses  a  steady  state  pitch  angle 
of  zero  radian  plus  or  minus  one-twentieth  radian.  In  this  configuration,  the  main  satellite 
bus  maintains  a  nadir  pointing  attitude.  The  region  of  stability  for  the  modulus  of  elasticity 
was  first  determined  while  holding  constant  the  dampening  coefficient  at  0.0007  Ns/m  for 
Satellite  I.  After  varying  the  modulus  of  elasticity  for  tether  one  between  zero  and  0.006 
N /mm2  and  the  modulus  of  elasticity  for  tether  two  between  zero  and  0.06  N /mm2,  the 
stable  range  was  found  as  depicted  in  Figure  5.21.  This  shows  that  stable  values  of  moduli 


Range  of  Moduli  of  Elasticity  for  Satellite  I 
(w/  average<0.05) 


Figure  5.21  Stability  Region  for  Moduli  of  Elasticity  for  Satl 
of  elasticity  exist  in  a  region  around  the  relationship  of  E2  equal  to  about  ten  times  E\. 
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This  relationship  may  be  explained  by  the  fact  that  the  distance  ”b”  is  a  little  over  ten  times 
the  value  of  ”a”.  Also,  the  main  satellite  bus  mass  is  about  ten  times  the  tip-mass.  This 
relationship,  however,  needs  more  support  before  a  correlation  can  be  confirmed.  However, 
there  are  a  few  data  points  shown  on  Figure  5.21  that  are  not  actually  stable  values  for 
moduli  of  elasticity.  These  are  quite  apparent  because  they  are  isolated  data  points  that  do 
not  follow  the  same  relationship  between  the  two  moduli  of  elasticity. 

5.4  Additional  Satellite  Results 

This  section  provides  the  results  of  Satellite  II  and  DumbSat.  The  specifications  for 
these  satellites  are  provided  in  Table  4.1.  These  results  serve  to  further  validate  the  program 
by  demonstrating  its  flexibility.  The  following  results  give  this  research  a  breadth  of  inputs 
upon  which  to  make  a  more  educated  conclusion.  Although  Satellites  I  and  II  are  very 
similar,  DumbSat  demonstrates  the  effectiveness  of  this  satellite  concept  and  simulation  for 
a  completely  different  satellite. 

5.4.1  Satellite  II  Results.  The  results  for  Satellite  II  are  summarized  in  this 
section.  Although  similar  analysis  has  been  performed  on  this  satellite’s  simulation  results, 
that  discussion  is  not  necessary  as  it  only  serves  to  back  up  the  conclusions  drawn  from 
Satellite  I.  For  Satellite  II,  the  region  of  stability  for  the  moduli  of  elasticity  is  shown  in 
Figure  5.22.  The  relationship  in  this  case  is  that  E2  is  equal  to  about  20  times  E\.  Also, 
the  distance  ”b”  is  about  24  times  the  value  of  ”a”,  and  the  main  satellite  bus  mass  is  about 
23  times  the  tip-mass.  Again,  there  is  not  enough  data  to  show  a  definite  correlation. 

For  Satellite  II,  the  lower  limit  for  stability  of  the  damping  coefficient  is  shown  in 
Figure  5.23.  For  the  nominal  values  of  moduli  of  elasticity  and  damping  coefficients  from 
Table  4.1,  the  resultant  graph  of  pitch  angle  is  shown  in  Figure  5.24.  For  this  graph  and  the 
others  in  this  section,  the  simulation  uses  1500  steps  and  an  interval  of  five  periods.  Figure 
5.25  is  the  graph  of  the  Sx  results,  and  Figure  5.26  is  the  graph  of  the  5y  results.  The  steady 
state  value  of  5x  is  about  28  meters  in  the  wake  of  the  satellite  while  oscillations  in  the  5y 
direction  are  less  than  a  meter  in  both  the  positive  and  negative  directions. 

As  shown  in  these  graphs,  the  results  for  Satellite  I  and  Satellite  II  are  very  similar. 
The  satellite  system  is  stable,  and  the  pitch  angle  dampens  to  approximately  zero  radian. 
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Range  of  Moduli  of  Elasticity  for  Satellite 
(w.  average<0.05) 
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Figure  5.22  Stability  Region  for  Moduli  of  Elasticity  for  Satll 


This  additional  case  supports  the  conclusions  drawn  from  Satellite  I.  These  results  demon¬ 
strate  and  support  the  conclusion  that  this  satellite  concept  is  capable  of  stabilizing  the 
attitude  to  a  relatively  constant,  nadir  orientation  within  a  limited  time  of  three  days. 
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Range  of  Damping  Coefficients  for  Satellite  II 
(w/E  1=0.0003  and  E2=0J00G) 
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Figure  5.23  Stability  Region  for  Damping  Coefficients  for  Satll 


0  vs.  Time  Step  for  Satll 


Figure  5.24  Nominal  Tether  Characteristics  for  Satll 
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5  x  vs.  Time  Step  for  Satll 


Figure  5.25  Nominal  Tether  Characteristics  for  Satll 


8  y  vs.  Time  Step  for  Satll 


Figure  5.26  Nominal  Tether  Characteristics  for  Satll 


5-22 


5.4-2  DumbSat  Results.  The  results  for  DumbSat  are  summarized  in  this  section. 
Although  similar  analysis  was  performed  on  this  satellite’s  simulation  results,  that  discussion 
is  not  necessary  as  it  only  served  to  back  up  the  conclusions  drawn  from  Satellite  I  and 
Satellite  II.  For  DumbSat,  the  region  of  stability  for  the  moduli  of  elasticity  is  shown  in 
Figure  5.27.  This  region  of  stability  shows  a  one-to-one  ratio  for  the  moduli  of  elasticity. 

Range  of  Moduli  of  Elasticity  for  DumbSat 
(for  average<0.05) 


Figure  5.27  Stability  Region  for  Moduli  of  Elasticity  for  DumbSat 

This  makes  sense  since  the  satellite  is  a  dumbbell  configuration  with  equal  end  masses  and 
equal  distances  to  each  end  mass  from  the  center  of  mass.  There  are  a  few  data  points 
shown  on  Figure  5.27  that  are  not  actually  stable  values  for  moduli  of  elasticity  because 
they  do  not  follow  the  same  relationship  between  the  two  moduli  of  elasticity. 

The  region  of  stability  for  the  damping  coefficients  is  shown  in  Figure  5.28.  For  the 
nominal  values  of  moduli  of  elasticity  and  damping  coefficients  from  Table  4.1,  the  resulting 
graph  of  pitch  angle  is  shown  in  Figure  5.29  For  this  output  and  the  following  graphs  in  this 
section,  the  simulation  uses  1500  steps  and  an  interval  of  five  periods.  Figure  5.31  shows 
the  graph  of  Sx  results,  and  Figure  5.31  is  the  graph  of  5y  results.  The  steady  state  value 
of  5x  is  about  28  meters  in  the  wake  of  the  satellite.  Oscillations  in  the  Sy  direction  are  less 
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0  (radians) 


Range  of  Damping  Coefficient  for  DumbSat 
(W  £^2=01)0317) 


ci  (Ns/hi) 


Figure  5.28  Stability  Region  for  Damping  Coefficients  for  DumbSat 


0  vs.  Time  Step  for  DumbSat 


Figure  5.29  Nominal  Tether  Characteristics  for  DumbSat 
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8  x  (meters) 


8  x  vs.  Time  Step  for  DumbSat 


Figure  5.30  Nominal  Tether  Characteristics  for  DumbSat 

8  y  vs.  Time  Step  for  DumbSat 


than  a  meter  in  both  the  positive  and  negative  directions. 
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This  final  case  supports  the  conclusions  drawn  from  Satellites  I  and  II.  These  results 
demonstrate  and  support  the  conclusion  that  this  satellite  concept  is  capable  of  stabilizing 
the  attitude  to  a  relatively  constant,  nadir  orientation  within  a  limited  time  of  three  days. 
This  case  also  shows  that  the  satellite  concept  can  be  used  for  different  gravity  gradient 
configurations,  largely  different  from  Satellites  I  and  II.  Although  stability  is  achieved,  this 
configuration  has  a  much  larger  settling  time  than  those  of  Satellites  I  and  II  and  so  would 
not  normally  be  chosen  over  Satellites  I  or  II. 
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VI.  Conclusions 


The  main  objective  of  this  research  was  to  determine  the  feasibility  of  using  a  gravity 
gradient  boom  and  a  tethered  balloon  for  passive  attitude  stabilization.  This  was  to  be 
accomplished  by  developing  and  testing  a  dynamical  model  for  the  satellite  system  through 
a  computer  program  that  shows  the  progress  of  the  pitch  angle  over  time.  A  literature 
review  was  conducted  to  understand  the  components  involved  in  this  satellite  concept.  The 
literature  review  also  discussed  past  research  that  used  similar  concepts  for  passive  attitude 
stabilization.  Equations  of  motion  were  calculated  for  the  three  mass  system  that  models 
the  satellite  concept  under  study.  The  program  was  validated  by  calculating  the  frequency 
of  oscillations  and  comparing  them  to  predicted  values  calculated  from  the  moments  of 
inertia.  The  simulation  program  was  run  for  three  different  satellites.  A  range  of  stability 
for  the  moduli  of  elasticity  and  damping  coefficients  was  found  for  each  satellite. 

The  simulations  run  show  that  a  gravity  gradient  satellite  with  a  tethered  balloon  can 
achieve  attitude  stabilization  over  time.  While  these  are  limited  cases,  this  basic  simulation 
demonstrates  the  feasibility  of  using  this  satellite  concept  for  aerostabilization.  The  gravity 
gradient  portion  of  the  attitude  control  system  provides  basic  stabilization,  oscillating  about 
its  equilibrium  point.  The  tethered  balloon  dampens  the  motion  by  dissipating  energy 
through  the  tethers. 

The  system  was  limited  to  a  circular  orbit  about  a  spherical  Earth  and  subjected 
to  only  drag  and  gravity  perturbations.  The  simulation  also  assumed  that  motion  was 
restricted  to  the  orbital  plane.  The  simulation  also  focused  on  modelling  only  microsatellites. 
Further  study  should  include  non-circular  orbits  and  include  more  orbital  perturbations, 
like  an  oblate  Earth.  This  would  make  it  necessary  to  include  motion  outside  of  the  orbital 
plane  and  determine  the  satellite  concept’s  effectiveness  in  damping  motion  in  the  roll  and 
yaw  directions.  Other  orbital  factors  could  be  analyzed  to  determine  its  impact  on  attitude 
stability  for  the  satellite  concept.  For  example,  the  maximum  altitude  with  effective  amounts 
of  aerodynamic  drag  should  be  determined. 

Once  the  simulation  is  updated  to  include  out-of-plane  motion,  the  model  should  be 
run  to  determine  the  stable  range  for  other  satellite  characteristics.  The  size  of  the  tethered 
balloon  should  also  be  varied  to  determine  the  minimum  size  that  would  create  enough  of  a 
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drag  force.  The  relationship  between  mission  lifetime  and  attitude  settling  time  should  be 
analyzed  to  determine  an  optimum  trade-off.  The  length  of  the  tethers  and  gravity  gradient 
boom  could  also  be  analyzed  to  determine  minimum  values  under  different  conditions.  The 
initial  conditions  should  also  be  varied  to  determine  if  the  satellite  concept  can  dampen 
out  unwanted  motion  despite  its  initial  orientation.  Since  there  are  so  many  variables  to 
include,  a  graphical  user  interface  should  be  created  that  can  make  it  easier  to  vary  certain 
quantities  and  provide  graphical  and  numerical  results  for  analysis. 

Further  research  should  also  include  analysis  and  inspection  of  possible  materials  for 
the  tethers  and  balloon.  A  mechanism  for  deployment  of  the  gravity  gradient  boom  and 
tethered  balloon  would  have  to  be  designed.  If  all  further  research  supports  the  feasibility 
and  effectiveness  of  this  satellite  concept  for  attitude  control,  a  test  flight  would  eventually 
have  to  be  conducted  to  verify  the  accuracy  of  the  simulation. 


6-2 


Appendix  A.  Thesis. m 


function  Thesis; 


0/ 0/ 0/ Of  0/  Of  0/0/ 0/0/ 0/0/0/  0/0/  0/0/  0/0/  0/0/  0/0/  0/0/  0/0/  0/0/  0/0/  0/0/  0/0/  0/0/  0/0/  0/0/  0/0/  0/0/  0/0/  0/0/  0/0/  0/0/  0/0/  0/0/  0/0/  0/0/  0/0/  0/0/  0/0/  0/0/  0/0/  0/0/0/ 

/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 


°/o°/o 

°/o°/o  PROCEDURE  Thesis 


11 

11 

11 

11 

11 

11 

11 

11 

11 

11 

11 

11 

11 

11 

11 

11 

11 

11 

11 

11 

11 

11 

11 

11 

11 

11 


This  procedure  uses  a  numerical  integrator  for  orbital  prediction  of 
a  three  body  system. 


Author  :  2Lt  Ernest  Maramba  AFIT/ENY  937-369-6956  27  Oct  2004 


Algorithm:  1.  Initialize  constants  and  variables 

2.  Declare  initial  conditions  and  satellite  characteristics 

3.  Call  RK4  and  RK4forRandV  which  predicts  generalized 
coordinates,  R  for  COM,  and  V  for  COM  at  next  time  step 

4.  Increment  time  step  and  loop  to  step  3  for  time  interval 

5.  Calculate  generalized  coordinates,  R  for  COM,  and  V  for 
COM  for  any  time  left  over 

6.  Plot  generalized  coordinates 


Locals 

Derivtype 

time 

nsteps 

leftover 

lamone 

lamtwo 

X 

Y 

steps 

period 


-  Defines  which  equations  of  motion  to  use 

-  Time  that  is  propagated  forward  in  Julian  time 

-  Number  of  steps  to  propagate  over 

-  Any  amount  of  time  not  integrated  over  by  initial  loop 

-  Length  of  tether  one  at  (time) 

-  Length  of  tether  two  at  (time) 

-  State  vector  for  R  and  V 

-  State  vector  for  generalized  coordinates 

-  Defines  number  of  steps  to  iterate  over 

-  Period  of  the  orbit 


A-l 


n 

interval 

Total  time  to  iterate  over 

n 

stepsize 

Stepsize  for  iterations 

n 

n 

Constants 

n 

Rad 

-  Conversion  for  degrees  to  radians 

n 

Deg 

-  Conversion  for  radians  to  degrees 

n 

TwoPI 

-  Pi  times  two 

n 

MU 

-  Constant 

n 

a 

-  Distance  from  main  satellite  bus  to  COM 

(m) 

n 

b 

-  Distance  from  tip-mass  to  COM  (m) 

n 

lonenot 

-  Taut  length  of  top  tether  (m) 

n 

ltwonot 

-  Taut  length  of  bottom  tether  (m) 

n 

ri 

-  Radius  of  circular  orbit  of  system  (km) 

n 

vi 

-  Velocity  of  circular  orbit  of  system  (km/sec) 

n 

mone 

-  Mass  of  tip-mass  (kg) 

n 

mtwo 

-  Mass  of  main  satellite  (kg) 

n 

mthree 

-  Mass  of  main  balloon  (kg) 

n 

JDstart 

-  Starting  Julian  date 

n 

JDstop 

-  End  Julian  date 

n 

n 

Coupling 

n 

Data 

-  Provides  global  constants 

n 

JulianDay 

-  Converts  year,  month,  etc.  into  a  Julian 

Date 

n 

RK4f orRandV 

-  Propagates  R  and  V  forward  in  time 

n 

RK4 

-  Propagates  delx,  dely,  and  theta  forward 

in  time 

n 

o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 
/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 

°/o°/o  Which  Satellite? 
question=l ; 

°/0°/o  Open  file  with  globals  and  declare  them 

Data (question)  global  Rad  Deg  MU  TwoPI  a  b  lonenot  ltwonot  ri  vi 
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JDstart  JDstop  mone  mtwo  mthree 


11  Initialize  Derivtype  -  2  means  the  com  is  only  subject  to  2-body 
11  motion 

Derivtype= Jynnnnnnnn2J ; 

11  Initialize  state  vector  Y  (del  x,  del  y,  theta) 

11  -  Y  is  in  terms  of  the  local  frame 

theta=pi/6 ;  delx=b*sin(theta) -lonenot*cos (theta- (pi/6) ) ; 

dely=b*cos (theta) +lonenot *sin (theta- (pi/6) ) ;  delxdot=0 ;  delydot=0 ; 

thetadot=0; 

Y  =  [delx;  dely;  theta;  delxdot ;  delydot;  thetadot] ; 

11  Initialize  state  vector  X  (satellite  COM) 

11  -  X  is  in  terms  of  the  inertial  GCE  frame 

11  -  Initialize  input  C0EJs  and  convert  units  11 

tempa=ri;  tempe=0.00;  tempi=0;  tempomega=0;  tempargp=0;  tempnu=0; 
tempu=0;  templ=0;  tempcappi=0; 

11  Convert  units  11 

tempi=tempi*Rad ;  tempomega=tempomega*Rad ;  tempargp=tempargp*Rad ; 
t empnu=t empnu*Rad ;  t empu=tempu*Rad ;  t empl=t empl*Rad ; 
tempcappi=tempcappi*Rad; 

11  Calculate  semi-latus  rectum  11 
tempp=tempa* (l-tempe~2) ; 

11  Calculate  R  and  V  11 
[R,V]  = 

RandV (tempp , tempe , tempi , tempomega,tempargp , tempnu,tempu, tempi ,tempcappi) 
X  =  [R(l : 3) ; V(1 : 3)]  ; 
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11  Define  time  increment  variable  and  number  of  steps 
steps=400;  period=TwoPI*sqrt (ri~3/MU) ;  interval=2*period; 
JDstop=JDstart+interval/86400 ;  stepsize=interval/steps ; 
time= JDstart ;  nsteps= ( JDstop- JDstart ) / (stepsize/86400) ; 


11  Initialize  El,  E2, 

cl,  and  c2 

f actor=l ; 

if  question==l 

11  Specs  for  Sat  I 

Eone=0. 002065; 

n 

Modulus 

Etwo=0.0202; 

n 

Modulus 

cone=0 . 0007*f actor ; 

n 

Damping 

ctwo=0 . 0007*f actor ; 

n 

Damping 

elseif  question==2 

11  Specs  for  Sat  II 

Eone=0 . 0003; 

n 

Modulus 

Etwo=0 . 006 ; 

n 

Modulus 

cone=0 . 0007 ; 

n 

Damping 

ctwo=0 . 0007 ; 

n 

Damping 

elseif  question==3 

11  Specs  for  DumbSat 

Eone=0 . 00317 ; 

n 

Modulus 

Etwo=0 . 00317 ; 

n 

Modulus 

cone=0.0016; 

u 

Damping 

ctwo=0.0016; 

n 

Damping 

end 

of  Elasticity  for  tether  one  (N/mnT2) 
of  Elasticity  for  tether  two  (N/mnT2) 
coefficient  for  tether  one  (N*s/m) 
coefficient  for  tether  two  (N*s/m) 

of  Elasticity  for  tether  one  (N/mnT2) 
of  Elasticity  for  tether  two  (N/mnT2) 
coefficient  for  tether  one  (N*s/m) 
coefficient  for  tether  two  (N*s/m) 

of  Elasticity  for  tether  one  (N/mnT2) 
of  Elasticity  for  tether  two  (N/mnT2) 
coefficient  for  tether  one  (N*s/m) 
coefficient  for  tether  two  (N*s/m) 


11  Propagate  X  and  Y  forward 
for  i=l:nsteps 

11  Save  Y  data  to  be  plotted  in  the  body  frame 
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xdata(i)=Y(l) ; 
ydata(i)=Y(2) ; 
thetadata(i)=Y(3) ; 

°/o°/o  Calculate  Y  at  (time)  by  calling  RK4 

[Y,  Rot]  =  RK4  (time,  stepsize,  Derivtype,  Y,  X,  Eone,  Etwo, 
cone,  ctwo) ; 

11  Calculate  X  at  (time)  by  calling  RK4forRandV 
X  =  RK4forRandV  (  time,  stepsize,  Derivtype,  X,  Y,  Rot  ); 

11  Increment  time 
time=time+(stepsize/86400) ; 
end  aver  =  sum(thetadata(l : f ix(nsteps) ) ) /f ix(nsteps) ; 
maxt=max(thetadata(f ix(nsteps)-100 : f ix(nsteps) ) ) ; 
mint=min(thetadata(f ix(nsteps)-100 : f ix(nsteps) ) ) ; 

11  Propagate  after  any  time  left  over  11 
lef tover= ( JDstop-time) *86400 ; 

[Y,  Rot]  =  RK4  (time,  leftover,  Derivtype,  Y,  X,  Eone,  Etwo, 
cone ,  ctwo) ; 

X  =  RK4forRandV  (  time,  leftover,  Derivtype,  X,  Y,  Rot); 

11  Plot  generalized  coordinates  in  local  coordinate  frame 

subplot (2,2, 1 :2) 

plot (thetadata) 

xlabel(’Time  Steps’) 

ylabel ( ’ \theta  (radians) ’ ) 

if  question==l 

title (’ \theta  vs.  Time  Steps  for  Sat  I’) 
elseif  question==2 
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title ( ’ \theta  vs.  Time  Steps  for  Sat  II’) 
elseif  question==3 

title ( ’ \theta  vs.  Time  Steps  for  DumbSat’) 

end 

subplot (2,2,3) 
plot (xdata) 
xlabel(’Time  Steps’) 
ylabel ( ’ \delta  x  (meters)’) 
if  question==l 

title (’ \delta  x  vs.  Time  Steps  for  Sat  I’) 
elseif  question==2 

title (’ \delta  x  vs.  Time  Steps  for  Sat  II’) 
elseif  question==3 

title (’ \delta  x  vs.  Time  Steps  for  DumbSat’) 

end 

subplot (2,2,4) 
plot (ydata) 
xlabel(’Time  Steps’) 
ylabel ( ’ \delta  y  (meters)’) 
if  question==l 

title (’ \delta  y  vs.  Time  Steps  for  Sat  I’) 
elseif  question==2 

title (’ \delta  y  vs.  Time  Steps  for  Sat  II’) 
elseif  question==3 

title (’ \delta  y  vs.  Time  Steps  for  DumbSat’) 

end 
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Appendix  B.  RK^.m 


function  [Y2,  Rot]  =  RK4  (time,  stepsize,  Derivtype,  Y,  X,  Eone, 
Etwo,  cone,  ctwo) ; 


0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/ 0/ 
/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 

°/o°/o 

°/o°/o  function  RK4 

°/o°/o 

°/0°/o  This  is  a  fourth  order  Runge-Kutta  integrator  for  a  6  dimension 
11  First  Order  differential  equation.  It  is  for  a  satellite 
11  equation  of  motion.  The  user  must  provide  an  external  function 
11  containing  the  system  Equations  of  Motion. 

11  The  integration  is  done  for  one  time  step  only. 

°/o°/o 

11  This  program  was  altered  in  order  to  be  used  for  this  thesis. 

11  Algorithm  :  Evaluate  each  term  depending  on  the  derivtype 
11  Find  the  final  answer 
11 

11  Author  :  Capt  Dave  Vallado  USAFA/DFAS  719-472-4109  5  Jun  1991 
11  In  Ada  :  Dr  Ron  Lisowski  USAFA/DFAS  719-472-4110  12  Jan  1996 
11  In  MatLab  :  Dr  Ron  Lisowski  USAFA/DFAS  719-333-4109  16  Jan  2002 
11  Altered  by  :  2Lt  Ernest  Maramba  AFIT/ENY  937-369-6956  8  Nov  2004 
11 

11  Inputs  : 

11  time  -  Initial  Time  Julian  Date  days  since  4713  B.C. 

11  stepsize  -  Step  size  sec 

11  DerivType  -  String  containing  YN  for  tension  forces  "YYNYNYNY2" 

11  Y  -  State  vector  of  delx,  dely,  theta  at  initial  time  m,  m/sec 
11  X  -  State  vector  of  R  and  V  at  initial  time  km,  km/sec 
11 

11  Outputs  : 
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11  Y  -  State  vector  of  delx,  dely,  theta  at  new  time  m,  m/sec 
11 

11  Locals  : 

11  YDot  -  Derivative  of  State  Vector 

11  K1,K2,K3  -  Storage  for  values  of  state  vector  at  different  times 
11  (The  standard  Runge-Kutta  K  constants) 

11  TEMP  -  Storage  for  state  vector 

11  TempDate  -  Temporary  date  storage  half  way  between  dt  days  since 
11  4713  B . C . 

11 

11  Constants  : 

11  TwoPI  -  2  times  pi 

11  w  -  Angular  rate  of  the  satellite’s  com 
11  JDstart  -  The  initial  time  for  integration 
11 

11  Coupling  : 

11  TetherDeriv  function  for  Derivatives  of  E.O.M. 

11 

11  References  : 

11  Mathews,  "Numerical  Methods"  pg.  423-427 
11  BMW  pg.  414-415 
11 

0/0/  0/  0/  0/  0/  0/0/  0/0/0/  0/  0/ 0/ 0/0/ 0/0/ Of  Of  Of  Of  0/0/  0/0/0/ Of  Of  0/0/  0/0/  0/0/  Of  Of  0/0/  0/0/  0/0/  Of  Of  0/0/  0/0/0/ Of  Of  Of  0/0/ 0/0/ Of  Of  Of  0/0/  0/0/  0/0/0/ 

/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 


°/0°/o  Declare  constants 

global  TwoPI  MU  JDstart  a  b  lonenot  ltwonot  mone  mtwo  mthree 

11  calculate  angular  rate  of  orbit  (rad/sec) 
w=sqrt (MU/mag (X ( 1 : 3) ) ~3) ; 

11  Calculate  rotation  angle  and  perform  revolution  check 
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rotang=w* (time- JDstart) *86400;  rotang  =  revcheck(rotang ,  TwoPI) ; 


11  Calculate  rotation  matrix 

Rot=  [-sin (rotang)  cos(rotang)  0 ; cos (rotang)  sin(rotang)  0 ; 0  0  -1]  ; 

11  Redefine  the  tether  lengths  according  to  Y  at  (time) 
lamone  =  sqrt ( (b*sin(Y(3) )  -  Y(l))~2  +  (b*cos(Y(3))  -  Y(2))~2); 
lamtwo  =  sqrt ( (a*sin(Y(3) )  +  Y(l))~2  +  (a*cos(Y(3))  +  Y(2))~2); 

11  Determine  whether  to  add  the  tension  forces  of  each  tether 

11  -  The  tether  length  must  be  greater  than  its  initial  length  to 

11  include  tension  forces  on  each  mass  to  which  it  is  connected. 

if  lamone/lonenot  >  1.000 

Derivtype (3) =JyJ ; 

else 

Derivtype(3)=JnJ ; 

end  if  lamtwo/ltwonot  >  1.000 

Derivtype (4)  = J  y J  ; 

else 

Derivtype (4) =JnJ ; 
end 


11  Local  VARIABLES 

y.y.y.y.y.y.y.y.y.y.y.y.y.#/.y.y.  Evaluate  ist  Taylor  Series  Term  y.y.y.y.y.y.y.y.y.y.y.y.y.y.y.y. 

[YDot , ax2body , ay2body , at2body , axdrag , aydrag , atdrag , Tether lspringx , 
Tether ldamperx , Tether lspringy , Tether ldampery , Tether lspringtheta , 
Tether ldampertheta , Tether2springx , Tether2damperx , Tether2springy , 
Tether2dampery , Tether2springtheta , Tether2dampertheta]  = 
TetherDeriv (Derivtype ,  Y,  lamone,  lamtwo,  X,  Rot,  Eone,  Etwo, 
cone ,  ctwo) ; 
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7o7o7o7o7o7o7o7o7o7o7o7o7o7o7o7o  Update  Julian  Date  for  a  half  Dt  ninUUnUl- 

TempDate  =  time  +  stepsize  *  0.5  /  86400.0; 


ninnnnnni  Evaluate  2nd  Taylor  Series  Term  UnnUnUUn 
K1Y  =  stepsize  *  YDot; 

Kl(2)  =  stepsize  *  ax2body; 

Kl(3)  =  stepsize  *  ay2body; 

Kl(4)  =  stepsize  *  at2body; 

Kl(5)  =  stepsize  *  axdrag; 

Kl(6)  =  stepsize  *  aydrag; 

Kl(7)  =  stepsize  *  atdrag; 

Kl(8)  =  stepsize  *  Tether lspringx; 

Kl(9)  =  stepsize  *  Tether ldamperx; 

Kl(10)  =  stepsize  *  Tether lspringy; 

Kl(ll)  =  stepsize  *  Tether ldampery; 

Kl(12)  =  stepsize  *  Tetherlspringtheta; 

Kl(13)  =  stepsize  *  Tether ldamperthet a; 

Kl(14)  =  stepsize  *  Tether2springx; 

Kl(15)  =  stepsize  *  Tether2damperx; 

Kl(16)  =  stepsize  *  Tether2springy ; 

Kl(17)  =  stepsize  *  Tether2dampery ; 

Kl(18)  =  stepsize  *  Tether2springtheta; 

Kl(19)  =  stepsize  *  Tether2dampertheta; 

Temp  =  Y  +  0.5  *  K1Y ; 

11  Redefine  the  tether  lengths  according  to  Y  at  (time) 
lamone  =  sqrt ( (b*sin(Temp(3) )  -  Temp(l))~2  +  (b*cos (Temp(3) ) . . . 

-  Temp(2))~2);  lamtwo  =  sqrt ( (a*sin(Temp (3) )  +  Temp(l))~2  +. . . 
(a*cos(Temp(3))  +  Temp(2))~2); 

[YDot , ax2body , ay2body , at2body , axdrag , aydrag , atdrag , Tether lspringx , 
Tether ldamperx , Tether lspringy , Tether ldampery , Tetherlspringtheta , 
Tether ldampertheta , Tether2springx , Tether2damperx , Tether2springy , 
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Tether2dampery , Tether2spr ingthet a , Tether2damperthet a]  = 
TetherDeriv(Derivtype ,  Temp,  lamone,  lamtwo,  X,  Rot,  Eone,  Etwo 
cone,  ctwo) ; 


mmy.my.mn  Evaluate  3rd  Taylor  Series  Term  UnnUnUUn 
K2Y  =  stepsize  *  YDot ; 

K2(2)  =  stepsize  *  ax2body; 

K2(3)  =  stepsize  *  ay2body; 

K2(4)  =  stepsize  *  at2body; 

K2(5)  =  stepsize  *  axdrag; 

K2(6)  =  stepsize  *  aydrag; 

K2(7)  =  stepsize  *  atdrag; 

K2(8)  =  stepsize  *  Tether lspringx; 

K2(9)  =  stepsize  *  Tether ldamperx; 

K2(10)  =  stepsize  *  Tether lspringy; 

K2(ll)  =  stepsize  *  Tether ldampery; 

K2(12)  =  stepsize  *  Tetherlspringtheta; 

K2(13)  =  stepsize  *  Tether ldamperthet a; 

K2(14)  =  stepsize  *  Tether2springx; 

K2(15)  =  stepsize  *  Tether2damperx; 

K2(16)  =  stepsize  *  Tether2springy ; 

K2(17)  =  stepsize  *  Tether2dampery ; 

K2(18)  =  stepsize  *  Tether2springtheta; 

K2(19)  =  stepsize  *  Tether2dampertheta; 

Temp  =  Y  +  0.5  *  K2Y ; 

11  Redefine  the  tether  lengths  according  to  Y  at  (time) 
lamone  =  sqrt ( (b*sin(Temp (3) )  -  Temp(l))~2  +  (b*cos (Temp(3) ) . . . 

-  Temp(2))~2);  lamtwo  =  sqrt ( (a*sin(Temp (3) )  +  Temp(l))"2  +. . . 
(a*cos (Temp(3) )  +  Temp(2) ) "2) ; 

[YDot , ax2body , ay2body , at2body , axdrag , aydrag , atdrag , Tether lspringx , 
Tether ldamperx , Tether lspringy , Tether ldampery , Tetherlspringtheta , 
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Tether ldamperthet a , Tether2spr ingx , Tether2damperx , Tether2spr ingy , 
Tether2dampery , Tether2springtheta , Tether2dampertheta]  = 
TetherDeriv(Derivtype ,  Temp,  lamone,  lamtwo,  X,  Rot,  Eone,  Etwo, 
cone,  ctwo) ; 


nnnnnnnn  Evaluate  4th  Taylor 
K3Y  =  stepsize  *  YDot; 

K3(2)  =  stepsize  *  ax2body; 

K3(3)  =  stepsize  *  ay2body; 

K3(4)  =  stepsize  *  at2body; 

K3(5)  =  stepsize  *  axdrag; 

K3(6)  =  stepsize  *  aydrag; 

K3(7)  =  stepsize  *  atdrag; 

K3(8)  =  stepsize  *  Tether lspringx; 

K3(9)  =  stepsize  *  Tether ldamperx; 

K3(10)  =  stepsize  *  Tetherlspringy ; 

K3(ll)  =  stepsize  *  Tether ldampery; 

K3(12)  =  stepsize  *  Tetherlspringtheta; 

K3(13)  =  stepsize  *  Tether ldamperthet a; 

K3(14)  =  stepsize  *  Tether2springx; 

K3(15)  =  stepsize  *  Tether2damperx; 

K3(16)  =  stepsize  *  Tether2springy ; 

K3(17)  =  stepsize  *  Tether2dampery ; 

K3(18)  =  stepsize  *  Tether2springtheta; 

K3(19)  =  stepsize  *  Tether2dampertheta; 

Temp  =  Y  +  K3Y; 
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nnnnnnnn  update  Julian  Date  for  a  full  Dt  nnnnnnn- 

TempDate  =  time  +  stepsize  /  86400.0; 

7o°/o  Redefine  the  tether  lengths  according  to  Y  at  (time) 
lamone  =  sqrt ( (b*sin(Temp (3) )  -  Temp(l))~2  +  (b*cos (Temp(3) ) . . . 
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-  Temp(2))~2);  lamtwo  =  sqrt ( (a*sin(Temp(3) )  +  Temp(l))~2  +. . . 
(a*cos(Temp(3))  +  Temp(2))~2); 

[YDot , ax2body , ay2body , at2body , axdrag , aydrag , atdrag , Tether lspringx , 
Tether ldamperx , Tether lspringy , Tether ldampery , Tether lspringtheta , 
Tether ldampertheta , Tether2springx , Tether2damperx , Tether2springy , 
Tether2dampery , Tether2springtheta , Tether2dampertheta]  = 
TetherDeriv(Derivtype ,  Temp,  lamone,  lamtwo,  X,  Rot,  Eone,  Etwo, 
cone ,  ctwo) ; 

nnnnn-  Update  the  State  vector,  perform  integration  TLTLTL- 
K4Y  =  stepsize  *  YDot; 

K4(2)  =  stepsize  *  ax2body; 

K4(3)  =  stepsize  *  ay2body; 

K4(4)  =  stepsize  *  at2body; 

K4(5)  =  stepsize  *  axdrag; 

K4(6)  =  stepsize  *  aydrag; 

K4(7)  =  stepsize  *  atdrag; 

K4(8)  =  stepsize  *  Tether lspringx; 

K4(9)  =  stepsize  *  Tether lspringy ; 

K4(10)  =  stepsize  *  Tetherlspringtheta; 

K4(ll)  =  stepsize  *  Tether ldamperx; 

K4(12)  =  stepsize  *  Tether ldampery ; 

K4(13)  =  stepsize  *  Tetherldampertheta; 

K4(14)  =  stepsize  *  Tether2springx; 

K4(15)  =  stepsize  *  Tether2springy ; 

K4(16)  =  stepsize  *  Tether2springtheta; 

K4(17)  =  stepsize  *  Tether2damperx; 

K4(18)  =  stepsize  *  Tether2dampery ; 

K4(19)  =  stepsize  *  Tether2dampertheta; 

Y  =  Y  +  (  K1Y  +  2.0  *  (  K2Y  +  K3Y  )  +  K4Y)/  6.0  ; 
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for  i=2 : 19 


acceleration (i)=(  Kl(i)  +  2.0  *  (  K2(i)  +  K3(i)  )  +  K4(i))  /  6.0 
end 


11  Perform  revolution  check  on  theta 

I  Y (3)  =  revcheck  (Y(3) ,  TwoPI) ; 

II  For  some  reason,  revcheck  doesn’t  work. 

11  -  something  to  do  with  theta  being  close  to  zero 
Y (3)  =  Y (3)  -  TwoPI  *  f ix(Y (3)  /  TwoPI); 

11  Define  output  Y  vector 
Y2  =  Y; 
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Appendix  C.  TetherDeriv.m 


function 

[YDot , ax2body , ay2body , at2body , axdrag , aydrag , atdrag , Tether lspringx , 
Tether ldamperx , Tether lspringy , Tether ldampery , Tether lspringtheta , 

Tether ldamperthet a , Tether2spr ingx , Tether2damperx , Tether2spr ingy , 
Tether2dampery , Tether2spr ingthet a , Tether2damperthet a] 

=  TetherDeriv(Derivtype ,  Y,  lamone,  lamtwo,  X,  Rot,  Eone,  Etwo, 
cone,  ctwo) ; 

0/0/0/ Of 0/ 0/ 0/0/ 0/0/ 0/ 0/ 0/ 0/ 0/0/ 0/0/0/ Of  Of  Of  0/0/  0/0/  0/0/ Of  0/0/ 0/0/ 0/0/ Of  Of  0/0/  0/0/  0/0/ Of  Of  0/0/ 0/0/0/ Of  Of  Of  0/0/  0/0/  0/0/  Of  Of  0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 

/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 
°/o°/o 

°/o°/o  PROCEDURE  TetherDeriv 

°/o°/o 

11  This  procedure  calculates  Ydot . 

°/o°/o 

11  Algorithm:  1.  Call  global  variables 

11  2.  Call  Deriv  to  calculate  generalized 

11  coordinates  without  generalized  forces 

11  3.  Calculate  tether  tension  forces  and  dissipative 

11  forces 

11  4.  Sum  accelerations  to  get  final  generalized  coordinates 

°/o°/o 

11  Author  :  2Lt  Ernest  Maramba  AFIT/ENY  937-369-6956  27  Oct  2004 

°/o°/o 

11  Inputs  : 

11  Derivtype  -  Defines  which  equation  of  motion  to  use 

11  Y  -  State  vector  of  delx,  dely,  theta  at  initial  time  m,  m/sec 

11  lamone  -  Length  of  tether  one  at  (time) 

11  lamtwo  -  Length  of  tether  two  at  (time) 

11  X  -  State  vector  of  R  and  V  km,  km/sec 

11  Rot  -  Rotation  matrix 

11 
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11  Outputs  : 

11  YDot  -  Derivative  of  Y  state  vector 
11 

11  Locals  : 

11  R  -  Position  vector  km 
11  V  -  Velocity  vector  km/s 

11  ax  -  Acceleration  in  x  direction  from  perturbations 

11  ay  -  Acceleration  in  x  direction  from  perturbations 

11  az  -  Acceleration  in  x  direction  from  perturbations 

11  r  -  Magnitude  of  R 

11  temp (2)  -  Variables  that  increase  coding  efficiency 

11  RSun  -  Vector  to  sun  from  the  satellite 

11  RtAsc  -  Place  holder  for  calling  Sun 

11  Decl  -  Place  holder  for  calling  Sun 

11  dotsunr  -  Dot  product  of  the  RSun  and  R 

11  rho  -  Density  of  the  atmosphere 

11  VRel  -  Relative  velocity  vector  of  satellite 

11 

11  Globals  : 

11  MU  -  Mu  of  the  Earth 

11  OmegaEarth  -  Earth’s  rotational  rate 

11  a  -  Length  of  boom  from  COM  to  main  satellite  bus  (m) 

11  b  -  Length  of  boom  from  COM  to  tip-mass  (m) 

11  mone  -  Mass  of  tip-mass  (kg) 

11  mtwo  -  Mass  of  main  satellite  bus  (kg) 

11  mthree  -  Mass  of  balloon  (kg) 

11  lonenot  -  Taut  length  of  top  tether  (m) 

11  ltwonot  -  Taut  length  of  bottom  tether  (m) 

11  E  -  Modulus  of  Elasticity  of  both  tethers  (N) 

11 

11  Coupling  : 
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11  Data  -  Provides  constants 

11  Deriv  -  Calculates  Ydot  without  tension  forces 

11 

0/0/0/  Of Of Of Of  Of 0/0/0/ 0/ 0/ Of Of  Of 0/0/0/ Of  Of  Of  Of  Of  Of  Of  Of  Of  Of  Of  Of Of  Of Of  Of 0/ 0/ 0/0/ 0/0/0/ Of  Of  Of  Of  Of 0/0/0/ Of  Of  Of  0/0/  0/0/  Of  Of  Of  0/0/ 0/0/ 0/0/0/ 

/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 

°/o°/o  Declare  global  variables  11 

global  MU  OmegaEarth  a  b  mone  mtwo  mthree  lonenot  ltwonot 


11  calculate  radius  to  com 
r=mag(X(l : 3) ) ; 

11  calculate  angular  rate  of  orbit  (rad/sec) 
w=sqrt (MU/mag (X ( 1 : 3) ) "3) ; 

11  initialize  accelerations  11 
ax=0;  ay=0;  atheta=0; 

Tetherlspringx  =  0; 

Tether 1 damper x  =  0; 

Tetherlspringy  =  0; 

Tether ldampery  =  0; 

Tetherlspringtheta  =  0; 

Tether ldampertheta  =  0; 

Tether2springx  =  0; 

Tether2damperx  =  0; 

Tether2springy  =  0; 

Tether2dampery  =  0; 

Tether2springtheta  =  0; 
Tether2dampertheta  =  0; 


11  Call  the  function  Deriv  to  calculate  EOM  with  tethers  slack  H 
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[YDot , ax2body ,ay2body ,at2body ,axdrag,aydrag,atdrag]  =  Deriv  (Y,  X, 
Rot)  ; 


11  Calculate  acceleration  with  taut  1-one  11 

if  Derivtype(3)==,yJ , 

Tetherlspringx  =  Eone* (1/lonenot  -  1/lamone) * (b*sin(Y(3) )-Y(l) ) ; 
Tether ldamperx  =  -cone*(Y(4)-w*Y(2)-b*Y(6)*cos(Y(3))+w*b*cos(Y(3))) ; 
Tetherlspringy  =  Eone* (1/lonenot  -  1/lamone) * (b*cos (Y(3) )-Y (2) ) ; 
Tether ldampery  =  -cone*(Y(5)+w*Y(l)+b*Y(6)*sin(Y(3))-w*b*sin(Y(3))) ; 
Tetherlspringtheta  =  -  Eone* (1/lonenot  -  1/lamone) *( (b*sin(Y(3) ).. . 
-Y(l) ) *b*cos (Y (3) )  -  (b*cos (Y (3) )-Y (2) ) *b*sin(Y (3) ) ) ; 

Tether ldampertheta  =  cone*(b*cos(Y(3))*(Y(4)-w*Y(2)-b*Y(6) . . . 

*cos (Y (3) )+w*b*cos (Y (3) ) )  -  b*sin(Y(3) ) * (Y(5)+w*Y(l)+b*Y(6) . . . 

*sin(Y (3) )-w*b*sin(Y (3) ) ) ) ; 

ax=ax  +  Tetherlspringx  +  Tether ldamperx; 

ay=ay  +  Tetherlspringy  +  Tether ldampery ; 

atheta=atheta  +  Tetherlspringtheta  +  Tether ldampertheta; 

end 

11  Calculate  acceleration  with  taut  1-two  11 
if  Derivtype(4)==JyJ , 

Tether2springx  =  Etwo* (1/ltwonot  -  1/lamtwo) * (-a*sin(Y(3) ) -Y(l) ) ; 
Tether2damperx  =  -ctwo* (Y (4) -w*Y (2) +a*Y (6) *cos (Y (3) ) -w*a*cos (Y (3) ) ) ; 
Tether2springy  =  Etwo* (1/ltwonot  -  1/lamtwo) * (-a*cos (Y (3) ) -Y (2) ) ; 
Tether2dampery  =  -ctwo*(Y(5)+w*Y(l)-a*Y(6)*sin(Y(3))+w*a*sin(Y(3))) ; 
Tether2springtheta  =  -  Etwo* (1/ltwonot  -  1/lamtwo) *( (-a*sin(Y(3) ).. . 
-Y(l) ) *-a*cos (Y (3) )  +  (-a*cos(Y(3) )-Y(2) )*a*sin(Y(3) )) ; 
Tether2dampertheta  =  ctwo*(-a*cos(Y(3))*(Y(4)-w*Y(2)+a*Y(6) . . . 
*cos(Y(3))-w*a*cos(Y(3)))  -  b*sin(Y(3) ) * (Y(5)+w*Y(l) -a*Y(6) . . . 
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*sin(Y (3) )+w*a*sin(Y (3) ) ) ) ; 

ax=ax  +  Tether2springx  +  Tether2damperx; 

ay=ay  +  Tether2springy  +  Tether2dampery ; 

atheta=atheta  +  Tether2springtheta  +  Tether2dampertheta 

end 

°/0  Define  Acceleration  Terms  of  YDot  TL 
YDot(4:6)  =  [YDot (4)  +  ax/mthree;  ... 

YDot (5)  +  ay/mthree;  ... 

YDot (6)  +  atheta/(mone*b"2  +  mtwo*a~2)] ; 
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Appendix  D.  Deriv. m 


function  [YDot , ax2body , ay2body , at2body , axdrag , aydrag , atdrag]  = 
Deriv  (Y,  X,  Rot) 


0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/ 0/ 
/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 

°/o°/o 

°/o°/o  function  DERIV 

°/o°/o 

11  This  function  calculates  the  derivative  of  the  state  vector  for 
11  use  with  the  Runge-Kutta  algorithm. 

°/o°/o 

11  Author  :  2Lt  Ernest  Maramba  AFIT/ENY  937-369-6956  27  Oct  2004 

°/o°/o 

11  Algorithm:  1.  Call  global  variables 

11  2.  Call  Drag  program  to  calculate  the  forces  of  drag  and 
11  direction  of  drag 

11  3.  Calculate  accelerations  due  to  drag  and  gravity 
11  4.  Sum  accelerations  to  get  YDot 

°/o°/o 

11  Inputs  : 

11  X  -  State  Vector  for  R  and  V  km,  km/sec 
11  Y  -  State  Vector  for  generalized  coordinates 
11  Rot  -  Rotation  matrix 
11 

11  Outputs  : 

11  YDot  -  Derivative  of  State  Vector 
11 

11  Locals  : 

11  Rrel  -  Position  vector  of  mass  three 

11  r  -  Distance  from  center  of  Earth  to  satellite  center  of  mass 
11  w  -  Orbital  spin  rate 


D-l 


11  Vrel  -  Relative  velocity 


11  rho  -  Density 

11  Fdrag  -  Drag  force  on  balloon 

11 

11 

11  Constants  : 

11  MU  -  Global  constant 
11  OmegaEarth  -  Spin  rate  of  Earth 

11  a  -  Length  from  center  of  mass  to  main  satellite  (mass  two) 
11  b  -  Length  from  center  of  mass  to  tip  mass  (mass  one) 

11  mone  -  Tip  mass 

11  mtwo  -  Main  satellite  mass 

11  bcthree  -  Ballistic  coefficient  of  tethered  balloon 
11  w  -  Angular  rate  of  the  satellite’s  com 
11 

11  Coupling  : 

11  Atmos  -  Gets  density  of  atmosphere 
11  cross  -  Calculates  cross  product 
11  mag  -  Calculates  magnitude  of  a  vector 
11 

11  References  : 

11  None. 

11 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 
/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 


°/0°/o  Global  Constants 


global  MU  OmegaEarth  a  b  mone  mtwo  bcthree  w  r 


iiiiiiiiiiiiiiiii  calculate  the  drag  terms  yoyoyoyoyoyonyo7oyoyoyoo/oo/oo/oo/oo/oo/oo/o 
[adragl,  Vhatl,  adrag2,  Vhat2,  adrag3,  Vhat3]  =  Drag  (Y,  X,  Rot); 
°/0  adragl=0; 

°/0  adrag2=0; 
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°/0  adrag3=0; 


11  calculate  radius  to  com 
rcom=mag(X(l : 3) ) ; 

11  calculate  angular  rate  of  orbit  (rad/sec) 
w=sqrt (MU/mag (X ( 1 : 3) ) "3) ; 

ax2body  =  2*w*Y(5)  +  3*Y(2) *Y(1) / (rcom*1000) *w~2 ;  ay2body  = 

-2*w*Y (4)  +  3*Y(2)*w~2  +  9/2*Y(2) ~2/ (rcom*1000) *w~2  + 

3/2*Y(l) "2/ (r com* 1000) *w~ 2;  at2body  = 

(-mone* (3*w~2*b~2*cos (Y (3) ) *sin (Y (3) ) +3/2*w~2*b~3*sin (Y (3) ) / (rcom*1000) ) 


mtwo*(3*w~2*a~2*cos(Y(3))*sin(Y(3))-3/2*w‘''2*a~3*sin(Y(3))/(rcom*1000)))  .  .  . 
/(mone*b~2  +  mtwo*a~2) ;  axdrag  =  adrag3*Vhat3 (1) ;  aydrag  = 
adrag3*Vhat3(2) ;  atdrag  =  (-mone* (adragl*Vhatl (1) *b*cos (Y(3) )  - 
adragl*Vhat 1 (2) *b*sin(Y (3) ) )  -  mtwo* (adrag2*Vhat2 (1) *-a*cos (Y (3) ) 

+  adrag2*Vhat2(2)*a*sin(Y(3))))/(mone*b~2  +  mtwo*a~2) ; 


y.y.y.y.y.y.y.y.y.y.y.y.y.y.y.y.y.#/.y.y.  velocity  Terms  in  m/s  in  body  frame  y.y.y.y.y.y.y.y.y.y.y.y.y.y.y. 

YDot=  [Y (4) ;  ... 

Y (5) ;  ... 

Y (6) ;  ... 

y.y.y.y.y.y.y.y.y.y.y.y.y.#/.y.y.y.y.  Acceleration  Terms  in  m/s "2  in  body  frame  1111111 

1  2*w*Y(5)  +  adrag*Vhat (1) ;  ... 

°/0  -2*w*Y  (4)  +  3*Y(2)*w~2  +  adrag*Vhat  (2)  ;  ... 

°/0  (-mone*(3*w~2*b~2*cos(Y(3))*sin(Y(3))+3/2*w~2*b~3*sin(Y(3))/.  .  . 
(rcom*1000))  -  mtwo* (3*w~2*a~2*cos (Y(3) ) *sin(Y(3) ) -3/2*w~2*a~3 . . . 
*sin(Y(3))/(rcom*1000)))/(mone*b~2  +  mtwo*a~2)] ; 


1  The  following  contains  the  equations  of  motion  with  greater  order 
1  of  accuracy: 
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ax2body  + 
ay2body  + 
at2body  + 


axdrag;  . . . 
aydrag;  . .  . 
at drag] ; 
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Appendix  E.  Drag.m 

function  [adragl,  Vhatl,  adrag2,  Vhat2,  adrag3,  Vhat3]  =  Drag  (Y,  X,  Rot) 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 
/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 

7.7. 

7.7.  function  Drag 

7.7. 

7.7.  This  function  calculates  drag  for  each  mass. 

7.7. 

7.7.  Author  :  2Lt  Ernest  Maramba  AFIT/ENY  937-369-6956  27  Oct  2004 

7.7. 


7.7. 

7.7. 

7.7. 

7.7. 

7.7. 

7.7. 

7.7. 

7.7. 

7.7. 

7.7. 

7.7. 

7.7. 

7.7. 


Algorithm:  1.  Call  global  variables 

2.  Calculate  r  and  w  for  COM 

3.  For  ml:  calculate  r  and  v  in  km  and  km/s  and  in  inertial  frame 
calculate  atmospheric  density 

calculate  magnitude  of  air  drag  deceleration  in  m/s~2 
calculate  unit  vector  for  direction  of  drag  term 

4.  Repeat  step  3  for  m2  and  m3 


Inputs  : 

X  -  State  Vector  for  R  and  V  km,  km/sec 
Y  -  State  Vector  for  generalized  coordinates 
Rot  -  Rotation  matrix 


7.7.  Outputs  : 

7.7.  YDot  -  Derivative  of  State  Vector 

7.7. 


7.7.  Locals  : 

7.7.  Rrel  -  Position  vector  of  mass  three 

7.7.  r  -  Distance  from  center  of  Earth  to  satellite  center  of  mass 

7.7.  w  -  Orbital  spin  rate 
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11  Vrel  -  Relative  velocity 


11  rho  -  Density 

11  Fdrag  -  Drag  force  on  balloon 

11 

11 

11  Constants  : 

11  MU  -  Global  constant 
11  OmegaEarth  -  Spin  rate  of  Earth 

11  a  -  Length  from  center  of  mass  to  main  satellite  (mass  two) 
11  b  -  Length  from  center  of  mass  to  tip  mass  (mass  one) 

11  mone  -  Tip  mass 

11  mtwo  -  Main  satellite  mass 

11  bcthree  -  Ballistic  coefficient  of  tethered  balloon 
11  w  -  Angular  rate  of  the  satellite’s  com 
11 

11  Coupling  : 

11  Atmos  -  Gets  density  of  atmosphere 
11  cross  -  Calculates  cross  product 
11  mag  -  Calculates  magnitude  of  a  vector 
11 

11  References  : 

11  None. 

11 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 
/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 


°/0°/o  Global  Constants 


global  MU  OmegaEarth  a  b  mone  mtwo  bcone  bctwo  bcthree 


11  calculate  radius  to  com 
r=mag(X(l : 3) ) ; 

11  calculate  angular  rate  of  orbit  (rad/sec) 
w=sqrt (MU/mag (X(l : 3) ) "3) ; 
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7o7oyoyoyoyo0/o0/o0/o7o0/o0/o0/o0/o0/o0/o°/o  Calculate  the  drag  term  for  ml  70707070707070707070707070707070707070% 
%  Calculate  r  in  km  and  in  inertial  frame  by  adding  R  of  COM  to 
%  rotated  delx  and  dely 

Tempr= [b*sin(Y(3) ) ;  b*cos(Y(3));  0];  Tempr  =  Rot*Tempr/1000 ; 

Rrel= [X(l)+Tempr (1) ;X(2)+Tempr (2) ;X(3)+Tempr (3)] ; 


%  Calculate  v  in  km/s  and  in  inertial  frame  by  adding  V  of  COM  to 
%  rotated  delxdot  and  delydot 

Tempv=[b*cos(Y(3))*(Y(6)-w) ;  -b*sin(Y(3) ) *Y(6)  +  w*b*sin(Y(3) ) ; 

0] ;  Tempv  =  Rot *Tempv/ 1000 ; 

Vrel= [X(4)+Tempv(l) ;X(5)+Tempv(2) ;X(6)+Tempv(3)] . . . 

-cross ( [0 ; 0 ; OmegaEarth]  ,Rrel); 


%  Calculate  atmospheric  density  in  kg/nT3 
[rho]  =  ATMOS (Rrel)*(10~3) ; 


%  Calculate  magnitude  of  air  drag  deceleration  in  m/s~2 
adragl=-0 . 5*rho* ( (mag (Vrel) *1000) "2) / (bcone) ; 

%  Calculate  unit  vector  for  direction  of  drag  term 
%  ~  Rotated  into  body  frame 
Vtemp=Vrel/mag(Vrel) ;  Vhatl=Rot*Vtemp ; 


7o7o7o7o7o7o7o7o%7o7o7o7o7o7o7o7o  Calculate  the  drag  term  for  m2  nnUnnUniniU 
%  Calculate  r  in  km  and  in  inertial  frame  by  adding  R  of  COM  to  rotated 
%  delx  and  dely 

Tempr=[-a*sin(Y(3)) ;  -a*cos (Y(3) ) ;  0];  Tempr  =  Rot*Tempr/1000; 

Rrel= [X(l)+Tempr (1) ;X(2)+Tempr (2) ;X(3)+Tempr (3)] ; 
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°/0  Calculate  v  in  km/s  and  in  inertial  frame  by  adding  V  of  COM  to 
°/0  rotated  delxdot  and  delydot 

Tempv=[-a*cos(Y(3))*(Y(6)-w) ;  a*sin(Y(3) ) *Y(6)  -  w*a*sin(Y(3) ) ; 

0] ;  Tempv  =  Rot *Tempv/ 1000 ; 

Vrel=  [X(4)+Tempv(l) ;X(5)+Tempv(2) ;X(6)+Tempv(3)] . . . 

-cross ( [0 ; 0 ; OmegaEarth] ,Rrel); 

°/0  Calculate  atmospheric  density  in  kg/nT3 
[rho]  =  ATMOS (Rrel)*(10~3) ; 

°/0  Calculate  magnitude  of  air  drag  deceleration  in  m/s~2 
adrag2=-0 . 5*rho* ( (mag (Vrel) *1000) "2) / (bctwo) ; 

°/0  Calculate  unit  vector  for  direction  of  drag  term 
°/0  ~  Rotated  into  body  frame 
Vtemp=Vrel/mag(Vrel) ;  Vhat2=Rot*Vtemp ; 


nnnnnnnni  calculate  the  drag  term  for  m3  nunnunnuun 

°/0  Calculate  r  in  km  and  in  inertial  frame  by  adding  R  of  COM  to  rotated 
°/0  delx  and  dely 

Tempr=[Y(l);  Y(2);  0];  Tempr  =  Rot*Tempr/1000 ; 

Rrel= [X ( 1 ) +Tempr ( 1 ) ; X(2)+Tempr (2) ;X(3)+Tempr (3)] ; 


°/0  Calculate  v  in  km/s  and  in  inertial  frame  by  adding  V  of  COM  to 
°/0  rotated  delxdot  and  delydot 

Tempv=[Y(4)  -  w*Y(2);  Y(5)  +  w*Y(l);  0];  Tempv  =  Rot*Tempv/1000 ; 
Vrel= [X(4)+Tempv(l) ;X(5)+Tempv(2) ;X(6)+Tempv(3)] . . . 

-cross ( [0 ; 0 ; OmegaEarth] ,Rrel); 
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°/0  Calculate  atmospheric  density  in  kg/nT3 
[rho]  =  ATMOS (Rrel)* (1CT3) ; 

°/0  Calculate  magnitude  of  air  drag  deceleration  in  m/s~2 
adrag3=-0 . 5*rho* ( (mag (Vrel) *1000) "2) / (bcthree) ; 

°/0  Calculate  unit  vector  for  direction  of  drag  term 
°/0  -  Rotated  into  body  frame 
Vtemp=Vrel/mag(Vrel) ;  Vhat3=Rot*Vtemp ; 
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Appendix  F.  Atmos. m 


function  [adragl,  Vhatl,  adrag2,  Vhat2,  adrag3,  Vhat3]  =  Drag  (Y,X,  Rot) 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 
/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 

7.7. 

7.7.  function  Drag 

7.7. 

7.7.  This  function  calculates  drag  for  each  mass. 

7.7. 

7.7.  Author  :  2Lt  Ernest  Maramba  AFIT/ENY  937-369-6956  27  Oct  2004 

7.7. 


7.7. 

7.7. 

7.7. 

7.7. 

7.7. 

7.7. 

7.7. 

7.7. 

7.7. 

7.7. 

7.7. 

7.7. 

7.7. 


Algorithm:  1.  Call  global  variables 

2.  Calculate  r  and  w  for  COM 

3.  For  ml:  calculate  r  and  v  in  km  and  km/s  and  in  inertial  frame 
calculate  atmospheric  density 

calculate  magnitude  of  air  drag  deceleration  in  m/s~2 
calculate  unit  vector  for  direction  of  drag  term 

4.  Repeat  step  3  for  m2  and  m3 


Inputs  : 

X  -  State  Vector  for  R  and  V  km,  km/sec 
Y  -  State  Vector  for  generalized  coordinates 
Rot  -  Rotation  matrix 


7.7.  Outputs  : 

7.7.  YDot  -  Derivative  of  State  Vector 

7.7. 


7.7.  Locals  : 

7.7.  Rrel  -  Position  vector  of  mass  three 

7.7.  r  -  Distance  from  center  of  Earth  to  satellite  center  of  mass 

7.7.  w  -  Orbital  spin  rate 
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11  Vrel  -  Relative  velocity 


11  rho  -  Density 

11  Fdrag  -  Drag  force  on  balloon 

11 

11 

11  Constants  : 

11  MU  -  Global  constant 
11  OmegaEarth  -  Spin  rate  of  Earth 

11  a  -  Length  from  center  of  mass  to  main  satellite  (mass  two) 
11  b  -  Length  from  center  of  mass  to  tip  mass  (mass  one) 

11  mone  -  Tip  mass 

11  mtwo  -  Main  satellite  mass 

11  bcthree  -  Ballistic  coefficient  of  tethered  balloon 
11  w  -  Angular  rate  of  the  satellite’s  com 
11 

11  Coupling  : 

11  Atmos  -  Gets  density  of  atmosphere 
11  cross  -  Calculates  cross  product 
11  mag  -  Calculates  magnitude  of  a  vector 
11 

11  References  : 

11  None. 

11 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 
/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 


11  Global  Constants 


global  MU  OmegaEarth  a  b  mone  mtwo  bcone  bctwo  bcthree 


11  calculate  radius  to  com 
r=mag(X(l : 3) ) ; 

11  calculate  angular  rate  of  orbit  (rad/sec) 
w=sqrt (MU/mag (X(l : 3) ) "3) ; 
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7o7oyoyoyoyo0/o0/o0/o7o0/o0/o0/o0/o0/o0/o°/o  Calculate  the  drag  term  for  ml  70707070707070707070707070707070707070% 
%  Calculate  r  in  km  and  in  inertial  frame  by  adding  R  of  COM  to 
%  rotated  delx  and  dely 

Tempr= [b*sin(Y(3) ) ;  b*cos(Y(3));  0];  Tempr  =  Rot*Tempr/1000 ; 

Rrel= [X(l)+Tempr (1) ;X(2)+Tempr (2) ;X(3)+Tempr (3)] ; 


%  Calculate  v  in  km/s  and  in  inertial  frame  by  adding  V  of  COM  to 
%  rotated  delxdot  and  delydot 

Tempv=[b*cos(Y(3))*(Y(6)-w) ;  -b*sin(Y(3) ) *Y(6)  +  w*b*sin(Y(3) ) ; 

0] ;  Tempv  =  Rot *Tempv/ 1000 ; 

Vrel= [X(4)+Tempv(l) ;X(5)+Tempv(2) ;X(6)+Tempv(3)] . . . 

-cross ( [0 ; 0 ; OmegaEarth]  ,Rrel); 


%  Calculate  atmospheric  density  in  kg/nT3 
[rho]  =  ATMOS (Rrel)*(10~3) ; 


%  Calculate  magnitude  of  air  drag  deceleration  in  m/s~2 
adragl=-0 . 5*rho* ( (mag (Vrel) *1000) "2) / (bcone) ; 

%  Calculate  unit  vector  for  direction  of  drag  term 
%  ~  Rotated  into  body  frame 
Vtemp=Vrel/mag(Vrel) ;  Vhatl=Rot*Vtemp ; 


7o7o7o7o7o7o7o7o%7o7o7o7o7o7o7o7o  Calculate  the  drag  term  for  m2  nUnnUnnUUn 
%  Calculate  r  in  km  and  in  inertial  frame  by  adding  R  of  COM  to 
%  rotated  delx  and  dely 

Tempr=[-a*sin(Y(3)) ;  -a*cos (Y(3) ) ;  0];  Tempr  =  Rot*Tempr/1000; 

Rrel= [X(l)+Tempr (1) ;X(2)+Tempr (2) ;X(3)+Tempr (3)] ; 
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°/0  Calculate  v  in  km/s  and  in  inertial  frame  by  adding  V  of  COM  to 
°/0  rotated  delxdot  and  delydot 

Tempv=[-a*cos(Y(3))*(Y(6)-w) ;  a*sin(Y(3) ) *Y(6)  -  w*a*sin(Y(3) ) ; 

0] ;  Tempv  =  Rot *Tempv/ 1000 ; 

Vrel=  [X(4)+Tempv(l) ;X(5)+Tempv(2) ;X(6)+Tempv(3)] . . . 

-cross ( [0 ; 0 ; OmegaEarth] ,Rrel); 

°/0  Calculate  atmospheric  density  in  kg/nT3 
[rho]  =  ATMOS (Rrel)*(10~3) ; 

°/0  Calculate  magnitude  of  air  drag  deceleration  in  m/s~2 
adrag2=-0 . 5*rho* ( (mag (Vrel) *1000) "2) / (bctwo) ; 

°/0  Calculate  unit  vector  for  direction  of  drag  term 
°/0  ~  Rotated  into  body  frame 
Vtemp=Vrel/mag(Vrel) ;  Vhat2=Rot*Vtemp ; 


nnnnnnnni  calculate  the  drag  term  for  m3  nunnunnuun 

°/0  Calculate  r  in  km  and  in  inertial  frame  by  adding  R  of  COM  to 
°/0  rotated  delx  and  dely 

Tempr=[Y(l);  Y(2);  0];  Tempr  =  Rot*Tempr/1000 ; 

Rrel= [X ( 1 ) +Tempr ( 1 ) ; X(2)+Tempr (2) ;X(3)+Tempr (3)] ; 


°/0  Calculate  v  in  km/s  and  in  inertial  frame  by  adding  V  of  COM  to 
°/0  rotated  delxdot  and  delydot 

Tempv=[Y(4)  -  w*Y(2);  Y(5)  +  w*Y(l);  0];  Tempv  =  Rot*Tempv/1000 ; 
Vrel= [X(4)+Tempv(l) ;X(5)+Tempv(2) ;X(6)+Tempv(3)] . . . 

-cross ( [0 ; 0 ; OmegaEarth] ,Rrel); 
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°/0  Calculate  atmospheric  density  in  kg/nT3 
[rho]  =  ATMOS (Rrel)* (1CT3) ; 

°/0  Calculate  magnitude  of  air  drag  deceleration  in  m/s~2 
adrag3=-0 . 5*rho* ( (mag (Vrel) *1000) "2) / (bcthree) ; 

°/0  Calculate  unit  vector  for  direction  of  drag  term 
°/0  -  Rotated  into  body  frame 
Vtemp=Vrel/mag(Vrel) ;  Vhat3=Rot*Vtemp ; 
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Appendix  G.  RK^forRandV.m 


function  [X2]  =  RK4forRandV  (  IDate,  stepsize,  DerivType,  X,  Y,  Rot) 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 

°/o°/o 

°/0°/o  function  RK4forRandV 

°/o°/o 

°/0°/o  This  function  is  a  fourth  order  Runge-Kutta  integrator  for  a 
11  6  dimension  First  Order  differential  equation.  The  intended 
11  use  is  for  a  satellite  equation  of  motion.  The  user  must  provide 
11  an  external  function  containing  the  system  Equations  of  Motion. 

11  Notice  Julian  date  is  included  since  some  applications  in  PDERIV 
11  may  need  this.  The  LAST  position  in  DerivType  is  a  flag  for 
11  two-body  motion.  Two-Body  motion  is  used  if  the  10th  element  is 
11  set  to  "2",  otherwise  the  Yes  and  No  values  determine  which 
11  perturbations  to  use. 

11  The  integration  is  done  for  one  time  step  only. 

11 

11  Algorithm  :  Evaluate  each  term  depending  on  the  derivtype 
11  Find  the  final  answer 

11  Notice  the  4th  k  must  be  mult  by  Dt  since  kl-k3  did  so  in  assignval 
11  Also,  the  4th  k  is  left  as  xdot  since  it  was  just  calculated 
11 

11  Author  :  Capt  Dave  Vallado  USAFA/DFAS  719-472-4109  5  Jun  1991 
11  In  Ada  :  Dr  Ron  Lisowski  USAFA/DFAS  719-472-4110  12  Jan  1996 
11  In  MatLab  :  Dr  Ron  Lisowski  USAFA/DFAS  719-333-4109  16  Jan  2002 
11 

11  Inputs  : 

11  IDate  -  Initial  Time  Julian  Date  days  since  4713  B.C. 

11  Dt  -  Step  size  sec 
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11  DerivType  -  String  containing  YN  for  incl  perts  "YYNYNYNY2" 

11  BC  -  Ballistic  Coefficient  kg/m2 

11  X  -  State  vector  at  initial  time  km,  km/sec 

11 

11  Outputs  : 

11  X  -  State  vector  at  new  time  km,  km/ sec 
11 

11  Locals  : 

11  XDot  -  Derivative  of  State  Vector 

11  K1,K2,K3  -  Storage  for  values  of  state  vector  at  different  times 
11  (The  standard  Runge-Kutta  K  constants) 

11  TEMP  -  Storage  for  state  vector 

11  TempDate  -  Temporary  date  storage  half  way  between  dt  days  since 
11  4713  B . C . 

11 

11  Constants  : 

11  None. 

11 

11  Coupling  : 

11  Deriv  function  for  Derivatives  of  E.O.M. 

11  PDeriv  function  for  Derivatives  of  E.O.M.  with  Perturbations 
11 

11  References  : 

11  Mathews,  "Numerical  Methods"  pg.  423-427 
11  BMW  pg.  414-415 
11 

o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 

global  BC 


U  Local  VARIABLES 

nnnnnnmx  Evaluate  1st  Taylor  Series  Term  yX/X/X/X/X/X/.*/. 


G-2 


if  DerivType (10 : 10)  ==  ’2’  , 

[XDot]  =  Derivf orRandVC  X  ); 

else 

[XDot]  =  PDerivf orRandVC  IDate ,X, DerivType ,BC,  Y,  Rot  ); 

end 


0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 
/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 


TempDate  =  IDate 


Update  Julian  Date  for  a  half  Dt 
+  stepsize  *  0.5  /  86400.0; 


0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/0/ 
/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 


0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o°/o  Evaluate  2nd  Taylor 

K1  =  stepsize  *  XDot; 

Temp  =  X  +  0.5  *  Kl; 


Series  Term 


0/ 0/ 0/ 0/ 0/ 0/0/ 0/0/ 0/0/0/ 0/ 0/ 
/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 


if  DerivType (10 : 10)  ==  ’2’  , 

[XDot]  =  Derivf orRandV(  Temp  ); 

else 

[XDot]  =  PDerivf orRandV(  IDate , Temp, DerivType ,BC,  Y,  Rot  ); 

end 


0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o°/o  Evaluate  3rd  Taylor  Series  Term  llllllllllUTL 
K2  =  stepsize  *  XDot; 

Temp  =  X  +  0.5  *  K2; 


if  DerivType (10 : 10)  ==  ’2’  , 

[XDot]  =  Derivf orRandV(  Temp  ); 

else 

[XDot]  =  PDerivf orRandV(  IDate, Temp , DerivType, BC,  Y,  Rot  ); 

end 


nnnnnnnn  Evaluate  4th  Taylor  Series  Term  lUlllinUlllll 
K3  =  stepsize  *  XDot; 
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Temp  =  X  +  K3; 


nnnnnnnn  update  Julian  Date  for  a  full  Dt  nnnnnnn- 

TempDate  =  IDate  +  stepsize  /  86400.0; 


if  DerivType (10 : 10)  ==  , 

[XDot]  =  Derivf orRandVC  Temp  ); 

else 

[XDot]  =  PDerivf orRandVC  IDate, Temp , DerivType, BC,  Y,  Rot  ); 

end 


V/XV/XVhVh-  Update  the  State  vector,  perform  integration  TLTLTL- 
X2  =  X  +  (  K1  +  2.0  *  (  K2  +  K3  )  +  stepsize  *  XDot)  /  6.0; 
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Appendix  H.  DerivforRandV.m 


function  [XDot]  =  Derivf orRandV  (  X  ) 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 
/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 

n 

7.7.  function  DERIV 

7.7. 

7.7.  This  function  calculates  the  derivative  of  the  two-body  state  vector 

7.7.  for  use  with  the  Runge-Kutta  algorithm.  Note  time  is  not  needed. 

7.7. 

7.7.  Algorithm  :  Find  the  answer 

7.7. 

7.7.  Author  :  Capt  Dave  Vallado  USAFA/DFAS  719-472-4109  28  Aug  1989 

7.7.  In  Ada  :  Dr  Ron  Lisowski  USAFA/DFAS  719-472-4110  5  Jan  1996 

7.7.  In  MatLab  :  Dr  Ron  Lisowski  USAFA/DFAS  719-333-4109  14  Nov  2001 

7.7. 

7.7.  Inputs  : 

7.7.  X  -  State  Vector  km,  km/sec 

7.7. 

7.7.  Outputs  : 

7.7.  XDot  -  Derivative  of  State  Vector  km/sec,  km/se2 

7.7. 

7.7.  Locals  : 

7.7.  RCubed  -  Cube  of  R 

7.7.  MU_R3  -  Mu  /  R  cubed 

7.7. 

7.7.  Constants  : 

7.7.  None. 

7.7. 

7.7.  Coupling  : 

7.7.  None. 
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n 

7.7.  References  : 

7.7.  None. 

7.7. 

/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 

7.7.  Global  Constants 
global  MU 


7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.  Build  the  XDot  Vector  7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7. 
R=  sqrt (X ( 1 ) ~2+X(2) ~2+X(3) "2) ; 

RCubed=  R*R*R; 

MU_R3  =  -MU/RCubed; 

7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.  Velocity  Terms  7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7. 
XDot  =  [X(4) ;  ... 

X(5) ;  ... 

X(6) ;  ... 

7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.  Acceleration  Terms  7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7. 
X(l)  *  MU_R3 ;  . . . 

X(2)  *  MU_R3 ;  . . . 

X(3)  *  MU_R3] ; 
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